mcflugen/cem-rafem-coupling.ipynb

## cem-rafem-coupling.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_coast(spacing, z):\n",
    "    import matplotlib.pyplot as plt\n",
    "    from matplotlib.colors import LinearSegmentedColormap\n",
    "    \n",
    "    land = plt.cm.terrain(np.linspace(.4, 1, 128))\n",
    "    ocean = plt.cm.ocean(np.linspace(.5, .8, 128))\n",
    "    colors = np.vstack((ocean, land))\n",
    "    m = LinearSegmentedColormap.from_list('land_ocean', colors)\n",
    "    \n",
    "    (y, x) = np.meshgrid(np.arange(z.shape[0]) * spacing[0],\n",
    "                         np.arange(z.shape[1]) * spacing[1], indexing='ij')\n",
    "    \n",
    "    plt.pcolormesh(y * 1e-3, x * 1e-3, z, cmap=m, vmin=-100, vmax=100)\n",
    "    \n",
    "    plt.gca().set_aspect(1.) \n",
    "    plt.colorbar(orientation='horizontal').ax.set_xlabel('Elevation (m)')\n",
    "    plt.xlabel('Along shore (km)')\n",
    "    plt.ylabel('Cross shore (km)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building a coupled system of BMI-enabled models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the two models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/huttone/git/csdms/coupling/cmt/framework/bmi_setup.py:28: UserWarning: missing units for sediment_flux_flag\n",
      "  warnings.warn('missing units for {name}'.format(name=name))\n",
      "/Users/huttone/_z/conda/lib/python2.7/site-packages/matplotlib/__init__.py:1318: UserWarning:  This call to matplotlib.use() has no effect\n",
      "because the backend has already been chosen;\n",
      "matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "  warnings.warn(_use_error_msg)\n"
     ]
    }
   ],
   "source": [
    "from cmt.components import Cem, Rafem\n",
    "cem = Cem()\n",
    "raf = Rafem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cem.setup('_run_cem', number_of_cols=100, number_of_rows=300, grid_spacing=5000.)\n",
    "raf.setup('_run_rafem', number_of_columns=100, number_of_rows=300, row_spacing=5., column_spacing=5.,\n",
    "          rate_of_sea_level_rise=0.015, channel_discharge=10000.)\n",
    "cem.initialize('_run_cem/cem.txt')\n",
    "raf.initialize('_run_rafem/input.yaml')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *Rafem* model provides `land_surface__elevation` as an output item and *CEM* uses it as an input item. Great!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'land_surface__elevation'}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set(raf.get_output_var_names()) & set(cem.get_input_var_names())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up initial conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin, we set the initial conditions for the *CEM* model. We first initialize it's elevations to be that of the *Rafem* model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "z = raf.get_value('land_surface__elevation')\n",
    "cem.set_value('land_surface__elevation', z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then set the wave climate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cem.set_value(\"sea_surface_water_wave__height\", 2.)\n",
    "cem.set_value(\"sea_surface_water_wave__period\", 7.)\n",
    "cem.set_value(\"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\", 0. * np.pi / 180.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/huttone/_z/conda/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAADLCAYAAACS0c8UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucXVV99/HP98wlF0IYxsRJJLFBBTXFKqBAIYRBKKSW\nEvu0Al5QlNpX66NivSBgbcPraSnoI1rb0lYFiyBoVEDEC4nIIGgBwQQSQiB5SNQEMkkYhnDJZWbO\n7/ljrX1mnzNn5uxJzm3I7/16zWv2XvuyfvvMnPM7a1/WkpnhnHPOZZFrdADOOecmDk8azjnnMvOk\n4ZxzLjNPGs455zLzpOGccy4zTxrOOecyq3nSkLRR0sOSVki6P5Z1Slou6XFJyyR1pNa/WNI6SWsl\nnVbr+JxzzmVXj5aGAd1mdqSZHRPLLgKWm9nhwB1xHknzgbOB+cAi4CpJ3hpyzrkmUa8PZJXMnwlc\nG6evBd4epxcDN5rZgJltBNYDx+Ccc64p1Kul8VNJD0j6YCzrMrPeON0LdMXpVwCbUttuAg6pQ4zO\nOecyaK1DHSeY2VOSZgLLJa1NLzQzkzRWXybez4lzzjWJmicNM3sq/t4m6WbC6aZeSbPMbIuk2cDW\nuPpmYG5q8zmxrKBCgnHOOTcKMyu9VDBuqmWHhZKmAi1m9pykA4BlwKXAqcDTZnaFpIuADjO7KF4I\nv4GQWA4Bfgq8xlJBSjIzIz80BEA+n09NDzE4MBCmh4bI5/OF6eR3YXl+qHhZaj+DAwNF+wzr5QvT\nRctT9QwO7Cmse80N3+K8c94xcp09e4riKo2juO7yx5Y17vxQvuzrko7pjl8/xMlvPIKhweF951py\n5FrC/1YuF3+3iNa2XKEsvTyZbm1Lb5cbsV0yn95nMl3Yd9Hy4f3d8P3/4byzFoxYJ5kfWX86rpai\nOsvVXW67ZL61vWXEazHa6/JP//wjPvt3f1K0P7WG5YNmDKbeb4Nm7Irz6WXp36MtH4z72JX8X5Xu\nO7WfXfl80TZf/Ov/x+vfcT7VNpQfrn8oFUtSXlRW8rkzlLcRZRtv/m/mLn5fYf2i/Zepa7T6R10n\nwz7G2mey/q6ffou2U84eeYxljpvU/kjK88X7tX88Z0Q91SCpKkmj1i2NLuBmSUld3zSzZZIeAJZK\nOh/YCJwFYGZrJC0F1hD+7z9k3g2vc841jZomDTPbALypTHkfobVRbpvLgMtqGZdz+6NWKLQ2nNtb\n/gxEjRz5hiMaHUImh87uqrxSE3jj/LmVV2oCCxce1ugQKpox/8hGh5DJ9Ne+sdEhZJI79PcbHUJd\nedKokYmSNF41e1ajQ8jkjfNf2egQMjlp4eGNDqGimb9/VKNDyOSg1404SdGUWl41Md7r1eJJwznn\nXGaeNJxzzmXmScM551xmnjScc85l5knDOedcZpme05DUBZxA6FDwRWA18ICZ5WsYm3POuSYzZtKQ\ndDJhrIuXAb8m9BE1mdCV+WskfQf4gpntqHWgzjnnGq9SS+NtwAfN7LelCyS1AWcApwHfrUFszjnn\nmsyYScPMPjXGsgHg5qpH5JxzrmllvaZxMPBeYF5qGzOzj9YoLuecc00o691TPwJ+D3gYeAB4MP5U\nJKlF0gpJP4jznZKWS3pc0jJJHal1L5a0TtJaSaeN71Ccc87VWtZebieZ2cf3so4LCF2dHxjnLwKW\nm9nnJH06zidjaZwNzCeOpSHpcL9DyznnmkfWlsb1kv5K0uzYUuiU1FlpI0lzCBfTvwYkg3+cCVwb\np68l3IkFsBi40cwGzGwjsJ4wGJNzzrkmkbWlsRv4HPAZIPnmb8CrKmz3ReBTwPRUWZeZ9cbpXsJA\nTRCeAbk3td4mQovDOedck8iaND5JGHZ1e9YdSzoD2GpmKyR1l1vHzKzCmN8+ap9zzjWRrEljHbBz\nnPs+HjhT0tsIDwROl3Qd0CtplpltkTSb8MAgwGYgPdLOnFg2wpIlS7A4NvLChQs5ccGCcYbmnHMv\nbT09PfT09FR9v1mTxovASkl3Ek5VQYVbbs3sEuASAEknAZ80s3MlfQ54H3BF/H1L3ORW4AZJVxJO\nSx0G3F9u30uWLCE/NARAPp8vTDvnRudDve5furu76e7uLsxfeumlVdlv1qRxS/xJTheJ8Z86Sta/\nHFgq6XxgI3AWgJmtkbSUcKfVIPAhM/PTU85ViY8R7qoha9JYbWYPpAsk/WnWSszsLuCuON0HnDrK\nepcBl2Xdr3POufrKesvtVyS9IZmR9E7gs7UJyTnnXLPK2tL4C+C7kt4FnEjoUuSPahaVc865ppQp\naZjZE7F1cQvwG+B0M3uxppE555xrOpXG01hVUtRJOKV1nyQzsz+oWWTOOeeaTqWWRuaL3c455176\nKiWN7Wb2/FgrSDrQzJ6rYkzOOeeaVKW7p74v6QuSFko6ICmU9GpJ50taBiyqbYjOOeeaRaWWxqmE\nXmr/Gjg+9mw7CDwG/BB4r5ltqW2IzjnnmkWl4V6NkBx+WJ9wnHPONbOsD/c555xznjScc85l50nD\nOedcZpmThqQTJb0/Ts+UdGiF9SdLuk/SSkmrJS2J5Z2Slkt6XNIySR2pbS6WtE7SWkmn7eUxOeec\nq5FMSSN+4F8IXByL2oHrx9rGzHYBJ5vZm4A3AYskHQtcBCw3s8OBO+I8kuYDZwPzCbfxXiXJW0LO\nOddEsn4o/xmwGHgBwMw2AwdW2ijVP1U70EYYU+NM4NpYfi3w9ji9GLjRzAbMbCOwHjgmY3zOOefq\nIGvS2G1m+WQm/aDfWCTlJK0EeoFlZnY/0GVmvXGVXqArTr8C2JTafBNhBD/n3Cg27dkz6rJvLTux\njpG4/UXWpPEdSf8FdEj6K8Jppa9V2sjM8vH01BzgWElHlCw3xh4B0Efucy5l4549bBoYYOPu3azf\ntYs57e1sLJM4rr99AeecdrcnDld1FbtGlyTg28DrgOeAw4HPmtnyrJWY2bNxfPHTgV5Js8xsi6TZ\nwNa42mZgbmqzObFshCVLlmD50PBZuHAhJy5YkDUUl9LSOsTQYEujwwBgeudkci1iR9+uhsXQMWNK\nQ+svZ9OePQzGUY8HgXnt7QyaFYZt3ZXPM6+9nf/88fGcd/o9APz37Qt4z+n3MGjExHFSYX+LTump\n7wG4hunp6aGnp6fq+1WlYbhj0lhlZkeMueLI7WYAg2bWL2kKcDthfPBu4Gkzu0LSRUCHmV0UL4Tf\nQLiOcQjwU+A1peOExy7ZyQ8NAZDP51PTQwwODITpoSHyMbEUlg+llueHipel9jM4MFC0z7BevjBd\ntDxVz+DAnsK6pXUU1tmzpyiu0jiK6y5/bFnjTseRfl0AhgZ3YtbKpCnh5c3lBMCOvl3kWnLkWlRU\nnmsRrW05OmZMKczv6NtFLqfCuq1t6e1yRdslci0q2mcuJzq7phY+rNPT6Tha23Ijtkvvc2T96bha\nirYt3UdrW46OmSFhdMyYUpS8OrumFva1o2932WNIji9dZzoOtYblg2aFBJDM70oSghkbd+8utBzm\ntLWxfvduZrW1jdi2NGkkZV/7yQkAnHf6PQzGsmSdZPq2O7qplaH88LENpY4zKS8qK/ncGcpb2bL0\n+kXzZeoarf5R18mwj7H2WS7uSsdNan8k5fni/do/njOinmqQhJmp8ppjq9jSMDOT9KCkY+I1iaxm\nA9dKaiGcBvu2mf1I0r3AUknnAxuBs2I9ayQtBdYQvlR9qDRh7G929D3NtIM6iuahOClN73wZ/du3\n0TFj5oj1RzM0uBOpHbM8u3eK/NAQQ4MhqYQPz91lt5veOYn+7TsLH46dXVNjXCO/nU/vnFR2P9MP\nnsTzzw6fTumYOYW+3hcLH7z923fW/Rv/9M7JhfpC0hw+tv5tO2ltD62xjhmTi5LHQS+fwnNPVyfO\nTXv2jEgY89rbC0kli6SlMZZFp/TwkxomDvfSl3W41+OA90j6DfEOKkI+GXUQJjNbBRxVpryP0BFi\nuW0uAy7LGNNLWv/2bUzvfBl9vVsKrYPpnS8rLE/Ktj/1JB0zZtLXuyUkkG1bi9YrNbDnRdrapxa2\nH1Hvtp10zJxS9MEO0Nk1Mpn09b4YWgqzptK/bWfhFFMuFz5YO7um0L99F9M7JwHhW/eOZ3YXksL0\nzsn0b9s5IoYdfbsqJq/n+0e/ADweScsmaRmkjy2JeTiuEE/HzMnkcuLZrTs56OVT2Llj32JJWhiD\nZoXTT6+ZNKmoVeJcs8iaNE6Pv5P/4n1u4rjR7Ximj44ZM8nn83TMmAmMPD2V6OyaRX5oiM6uWQwO\nDBRaHtM7O4vW2/n8DnItoq19auX6+3YXfavu7JpCX+/OotNMaf3bdo748M3lcvT17hyRbJJTP8np\noPSHcuk+k+TV2RVOiT3fv4eOmZNjTFPLJpzxKG3lZLWjb3ch7uee3sWUzvgB/2L5RJxIX5+Y095e\nKJvnCcJNIJneLfG5iQ7CMxZ/ChwUy5wbU//22pxmShLVRDKnvb2QLJybqLI+EX4B4QnwmYTnKq6X\n9NFaBra/2tH3NNMP7qy84hjC9Y2+wvyz27cxZdp0Jk2p+DzmcBzP7C6cWurrrfyNPjmdk1WWVkLS\n4unr3cmOvt3hmsq24STU1/viuBNHcqqsGi2VxMBzA9jzg9jU0e9E27h7uLU1p72d9bt3h+sXnkTc\nBJP19NRfAsea2QsAki4H7gW+XKvAXor6tjwFDN/xVHrtIbkuMdr1hvGY3tlJ/7ZtDA0OcNCMmWVP\nbWUx2nWFekm3VMrFkiSO/u3DCSC5w+v5Z8P60zsnFd0BFRLRi6OeGttbenGI3VNyTNpTfKqp3Cmo\n5JrFeC50O9cMsiYNgPwo065E//atdMx4eWF+R9/TDA4M0DlrdtEtt/3btjKtoyNus61wXaJakqSU\n3Nb7UtW/ffiayvTOyYU7vDq7phYSRfoW2f7tIy98V8ukPcazbXBwfIdszg96a8K9pGRNGl8H7pN0\nE+Ei+NuBa2oW1QTXMePl9G8PzywmLYokUaSFi9ZbGRwYKFzwdnsnaXGkb9VNWh+tbfV9gPGgAXim\nLTxXcYjG873MueaX6T/azK6UdBewgHAH1XlmtqKmkVWw45m+sqd4Gq1/W7hVtmPGy0c8uFdOoTVQ\nhVNS+7tyt842ysF5AWqq+wwHK6/iXEXj+Rq0EtgStzFJrzSz39YmrLH1b99WuFjcv31b2QfadvQ9\nTT6fH5FUQnltks2Ovj46Zs4sPHznXDNpxROH23eZkoakjwD/QOgnKv2V+A21CKqSjhkzix5469+2\ntXAtYFpHR+EBt3w+X5ju691SWB9InT7KFz3TkKyXtGJ29PWVXV6aqMJF7H2768k555pd1pbGx4DX\nmtnTtQxmb6Uv+Pb1bqFj5ssLSaVj5svp37a18BBcklzSp4WSp6/TF6OHyztTiSRfWN6/bWtRf1Kh\nvDpPKTvnXLPKmjR+C+yoZSDVUu60U6VTUR0zZhb6byotz+eHYsLJF92FlE465S5yO9ds/NSUq4Yx\nk4akT8TJJ4AeSbcByddpM7MraxlcPVXr+QjnmpVf03DVUOmJ8AOBaYSWxnLCsK3TYnnFx4slzZV0\np6RHJK1OniKX1ClpuaTHJS2T1JHa5mJJ6yStlXTa3h6Yc8656huzpWFmS0rLYlfn08zs2Qz7HwD+\n1sxWSpoGPChpOfB+YLmZfU7Sp4GLgGRMjbOB+cQxNSQdnh5q1jnnXONk7XvqBknT49jgq4A1ki6s\ntJ2ZbTGzlXH6eeBRQjI4E7g2rnYt4WFBgMXAjWY2EDtEXE8YlMk551wTyNon9O+b2Q7Ch/uPgXnA\nueOpSNI84EjgPqDLzHrjol5CJ4gArwA2pTbbREgyzjnnmkDWpNEqqY2QNH5gZgMMj61RUTw19T3g\nAjN7Lr0sjs431r68RzfnnGsSWW+5/S/C0KwPAz+PrYYs1zSIyeZ7wHVmdkss7pU0y8y2SJpNeGgQ\nYDMwN7X5nFhWZMmSJVi8zXXhwoWcuGBBxsNwzrn9Q09PDz09PVXfb9a+p75Mqhv0OOzryZW2kyTg\namCNmX0ptehW4H3AFfH3LanyGyRdSTgtdRgwYlzyJUuWFG6PTboZd86NzW+33b90d3fT3d1dmL/0\n0kurst+96oIznlLK8j94AvAe4GFJSQeHFwOXA0slnU9owZwV97tG0lJgTdz/h2Jdzrl95M9puGqo\nab/NZnYPo183OXWUbS4DLqtZUM455/Za1gvhzrkJzlsZrhqyPqdxlqTpcfqzkm6WdFRtQ3POVZMP\nB+WqIWtL47NmtkPSAuAUwsXt/6hdWM4555pR1qSR3J50BvBVM7uN0A+Vc865/UjWpLFZ0lcI/UL9\nUNLkcWzrnHPuJSLrB/9ZwO3AaWbWDxwMfKpmUTnnnGtKWa+NzQJ+aGa7JJ0M/AHDHQ4655zbT2Rt\nadwEDEp6DaFLkTnADTWLyjnnXFPKmjTyZjYI/C/gX83sU8Ds2oXlnHOuGWVNGnskvQt4L3BbLGur\nTUjOOeeaVdak8QHgD4F/MrMNkl4FXFe7sJxzzjWjTEnDzB4BPgmslnQE8Dszu6KmkTnnnGs6WbsR\n6QYeB/4duApYJ+mkDNtdI6lX0qpUWaek5ZIel7RMUkdq2cWS1klaK+m0cR+Nc865msp6eupKwjMa\nC81sIXAa8MUM230dWFRSdhGw3MwOB+6I80iaT3h4cH7c5ipJ/gChc841kczDvZrZY8mMmT1Ohmc8\nzOxu4JmS4jMZfsbjWsIQsgCLgRvNbMDMNgLrgWMyxuecc64Osj7c96CkrwHXAwLeDTywl3V2mVlv\nnO4FuuL0K4B7U+ttIoze55xzrklkTRp/DXwY+Gicv5twbWOfmJlJGmtkPh+1zznnmkjFpCGpFXjI\nzF4HfKEKdfZKmmVmWyTNBrbG8s3A3NR6c2LZCEuWLMHyeQAWLlzIiQsWVCEs517afBCm/UtPTw89\nPT1V32+W6xKDkh6T9Htm9psq1Hkr8D7givj7llT5DZKuJJyWOgy4v9wOlixZQn4o9Naez+cL0865\n0fkY4fuX7u5uuru7C/OXXnppVfab9fRUJ/CIpPuBF2KZmdmZY20k6UbgJGCGpN8Bfw9cDiyVdD6w\nkdCDLma2RtJSYA3hf/tDZuanp5xzrolkTRqfLZnP9GFuZu8cZdGpo6x/GXBZxpicc87V2ZhJQ9Jh\nhLudekrKFwBP1TAu55xzTajScxpfAnaUKd8RlznnJgi/nuGqoVLS6DKzh0sLY9mhtQnJOVcLWc9F\nOzeWSkmjY4xlk6sZiHPOueZXKWk8IOmvSgslfRB4sDYhOeeca1aVWqwfA26W9G6Gk8TRwCTgz2oZ\nmHPOueYzZtKIT20fD5wMHEG41fY2M/tZPYJzzlWPXwh31ZDliXADfhZ/nHMTlD8R7qrBx6twzjmX\nmScN5/YT3spw1eBJw7n9hD+n4arBk4ZzzrnMmi5pSFokaa2kdZI+3eh4nHPODWuqpCGpBfg3YBEw\nH3inpNc3Nqq9s2LV6kaHkMkTT21pdAiZPLTmt40OIZO7fv54o0OoaNsjv250CJk8u3Zlo0PIZOiJ\nifFer5amShrAMcB6M9toZgPAt4DFDY5pr0yUpLHhqd7KKzWBh9b8rtEhZPLzn69rdAgVbV+zotEh\nZLLjsYcaHUIm+Q2PNDqEumq2pHEIkP502BTLnHPONYFmSxo+Up9zzjUxNdOIqpKOA5aY2aI4fzGQ\nN7MrUus0T8DOOTeBmJn2dR/NljRagceAU4AngfuBd5rZow0NzDnnHNBkz/uY2aCkDwO3Ay3A1Z4w\nnHOueTRVS8M551xza7YL4WNqlgf/JM2VdKekRyStlvTRWN4pabmkxyUtk9SR2ubiGPdaSafVOd4W\nSSsk/aBZ45TUIem7kh6VtEbSsU0a59/Gv/kqSTdImtQMcUq6RlKvpFWpsnHHJenoeGzrJP1LHWL8\nfPybPyTpJkkHNTLG0eJMLfuEpLykzmaNU9JH4mu6WlL6enB14jSzCfFDOF21HpgHtAErgdc3KJZZ\nwJvi9DTCdZjXA58DLozlnwYuj9PzY7xtMf71QK6O8X4c+CZwa5xvujiBa4EPxOlW4KBmi5Nw+/cT\nwKQ4/23gfc0QJ3AicCSwKlU2nriSsw73A8fE6R8Bi2oc4x8lrwlweaNjHC3OWD4X+AmwAehsxjgJ\nYx8tB9ri/MxqxzmRWhpN8+CfmW0xs5Vx+nngUcIHypmEDz/i77fH6cXAjWY2YGYbCX+wY+oRq6Q5\nwNuArwHJnRNNFWf8dnmimV0D4dqWmT3bbHFGrcDUeNPGVMINGw2P08zuBp4pKR5PXMdKmg0caGb3\nx/W+kdqmJjGa2XIzy8fZ+4A5jYxxtDijK4ELS8qaLc6/Af45fkZiZtuqHedEShpN+eCfpHmEbH8f\n0GVmySPWvUBXnH4FId5EPWP/IvApIJ8qa7Y4DwW2Sfq6pF9L+qqkA5otTjPbDHwB+C0hWfSb2fJm\nizNlvHGVlm+mvvF+gPBNlzKxNDRGSYuBTWb2cMmipooTOAxYKOleST2S3lztOCdS0mi6K/aSpgHf\nAy4ws+fSyyy09caKuebHI+kMYKuZrWC4lVEcRBPESfj2fhRwlZkdBbwAXFQURBPEKelgwrf3eYQ3\n2zRJ7ykKogniLFtp5bgaStJngD1mdkOjYyklaSpwCfAP6eIGhVNJK3CwmR1H+LK4tNoVTKSksZlw\nTjExl+IMWVeS2ggJ4zozuyUW90qaFZfPBrbG8tLY58SyWjseOFPSBuBG4K2SrmvCODcRvsX9Ks5/\nl5BEtjRZnKcCG8zsaTMbBG4C/rAJ40yM5++8KZbPKSmvebySziOcQn13qriZYnw14YvCQ/G9NAd4\nUFJXk8VJrPsmgPh+ykuaUc04J1LSeAA4TNI8Se3A2cCtjQhEkoCrgTVm9qXUolsJF0aJv29JlZ8j\nqV3SoYQm5P3UmJldYmZzzexQ4BzgZ2Z2bhPGuQX4naTDY9GpwCPAD5opTuA3wHGSpsT/gVOBNU0Y\nZ2Jcf+f4d9ihcOeagHNT29SEpEWEb8SLzWxXSexNEaOZrTKzLjM7NL6XNgFHxVN/TRNndAvwVoD4\nfmo3s+1VjbOaV/Nr/QP8MeFOpfXAxQ2MYwHhGsFKYEX8WQR0Aj8FHgeWAR2pbS6Jca8FTm9AzCcx\nfPdU08UJvBH4FfAQ4ZvSQU0a5xLCjQ+rCBeX25ohTkJL8klgD+Ha3/v3Ji7g6Hhs64Ev1zjGDwDr\nCMk4eR9d1cgYS+LcnbyWJcufIN491Wxxxv/H62K9DwLd1Y7TH+5zzjmX2UQ6PeWcc67BPGk455zL\nzJOGc865zDxpOOecy8yThnPOucw8aTjnnMvMk4abMCS9PXZL/dpU2bxyXVjXIZZuxa7m61DXlyQt\niNMb091yj3M/fyLp0upG5/Y3njTcRPJO4O74e0KTlOm9J+llwLFmdk8s2usHq8zsh8CfSpqyt/tw\nzpOGmxBi55AnAH9J6BKl3DqTY0+5D8fecrtj+XkKA/z8WGFAovTANOdLekzSfbF33X8ts9+TFAax\nWhH3Oy0umibpO3HAm+tT658S13tY0tWx25uklXC5pAeBd0g6TdIvJT0oaWns2bfUnwM/LhPTlHg8\n50v6PYWBdb4ej+V6SadKuice71tSm/YAZ4z5Yjs3Bk8abqJYDPzYzNYBT0s6qsw6/xsYMrM/ILRG\nrpU0KS57I3AW8AbgbEmHSHoF8HfAsYSE9FrKf5P/BPAhMzuS0IXMzlh+JHABYYCbV0k6XtJk4OvA\nWTGOVsIYB8R9bzezo4E7gM8Ap8T5BwmDZZU6Pi5LO5DQl9A3zexqQo+rrwb+L/C6+HOOmS0APkno\nPiLxAGHwHuf2iicNN1G8kzBSHvF3uVNUJwDXA5jZY4Q+jQ4nfFjfYWbPmdluQieD8wgDIt1lZv0W\neq39DuW7vP4F8EVJHyF0Oz0Uy+83syct9MWzkjAuyGsJPeGuj+tcCyxM7Ss5huMIyeaXklYA7wVe\nWabu2cC21LyA7wPXmNn1qfINZvZIjOURQlICWB2PNbGN0K27c3ultdEBOFdJvPB7MnCEJCMM/WuE\n3lFHrD7KbnanpocI//ulrYrRxhy5QtJtwJ8Av5B0+jj3mS57ITW93MzeNUq8iZ1A+hqEAfcQOu+8\nMVWejiVP6BQwmU6/zycz3FJybty8peEmgr8AvmFm8yx0T/1KYIOk0tMsdxPHZIjdQr+S0KNnuWRg\nhF51T5LUoTB8659T5vSUpFfHb/Gfi9uMdhrLCL0wz5P06lh2LnBXmXXvA05I1pN0gKTDyqz3KPCa\nkrK/B56R9O9l1q/kcEKPps7tFU8abiI4B7i5pOx7sTw9It1VQE7Sw4Qx5N9nYazksqPWmdmTwGWE\nMS7uATYAO8rUf4GkVZIeInyDTy5Ml9vnbkIX1d+JcQwC/1m6voWxm88Dboz7/SUhGZX6IdCdriJu\nfwEwRdLloxyfjTLdHffp3F7xrtHdfk3SAWb2Qmxp3ARcbWbfb3RcaZLuBs4ws2f3cT9dhIvnp1Yn\nMrc/8qTh9muSPk8YgW8ycLuZfazBIY0g6Rhgp5nt02klSW8mjMP9cHUic/sjTxrOOecy82sazjnn\nMvOk4ZxzLjNPGs455zKbcA/3xYe7nHPOjZOZjfbwa2YTLmkAcPISUOrYJciVzMNwWfI7vU1OI+fL\nbZtlm2S90m3LxTbaNqX7L91mrHrLHfto5cn2Y71eZWJokWgpmS/6nSueT8oqzZfbdq+2qUasZbYZ\nM4aSeo4+7vbCG6pVYnIuR2vcrjWWJfOTc7kRZa0Sk0vmi7YpmW+VsME8+byRHwrfpfJDRj5vDA7k\ni+YLy/PG4J4h8nkbsbywTZwv2mZgqGi+aP3UPpL9Dg7k4z7y5bcpWT9ZFsryI8rKbTO8fupYhvK0\ntLaRa2kBIJfLkWtpKcy3trWRy7WQa8nF5WFZa1tbmG9pKWxTKIvr5HK5wjqtbek6ipe3trcP15uL\n67TkaG1rL9rHyDpyqfWH6yjaJlcSa7m4UvtI4gBQ+jNgH/jpKeecc5l50nDOOZeZJw3nnHOZedJw\nzjmXmSfD/Ws0AAAJQElEQVQN55xzmXnScM45l5knDeecc5l50nDOOZeZJw3nnHOZedJwzjmXmScN\n55xzmXnS2Ftb1zc6Ati0ttERADD4xOpGh0D/2pWNDgGAtb94otEh8PO71zU6BB5+9HeNDgGAJ57c\n0ugQWLGq8e+PavKksbc8aRQMbWj8m+LZJkkaj/1iQ6ND4O57PGkknnjKk0a1edJwzjmXmScN55xz\nmclsYo1p5IMwOefc3qnGIEwTLmk455xrHD895ZxzLjNPGs455zJr2qQh6R2SHpE0JOmokmUXS1on\naa2k01LlR0taFZf9Sw1iepOkeyWtkPQrSW+pFFMtSPqIpEclrZZ0RSNiSNX5CUl5SZ31jkPS5+Pr\n8JCkmyQdVO8YUvUtinWtk/TpWtcX65wr6c74Plkt6aOxvFPSckmPS1omqaNO8bTE98YPGhGHpA5J\n343/E2skHduI10LS38a/xypJN0iaVOs4JF0jqVfSqlTZqHXu0/vDzJryB3gdcDhwJ3BUqnw+sBJo\nA+YB6xm+NnM/cEyc/hGwqMoxLQNOj9N/DNw5Rky5Gr0uJwPLgbY4P7PeMaRimQv8BNgAdDbgtfij\nZN/A5cDljXgtgJZYx7xY50rg9bV87WO9s4A3xelpwGPA64HPARfG8k8nr0sd4vk48E3g1jhf1ziA\na4EPxOlW4KAGxHAI8AQwKc5/G3hfreMATgSOBFalysrWua/vj6ZtaZjZWjN7vMyixcCNZjZgZhsJ\nB3yspNnAgWZ2f1zvG8DbqxxWnvCPCNABbB4jpmOqXHfib4B/NrMBADPb1oAYElcCF5aU1S0OM1tu\nZvk4ex8wp94xRMcA681sY/y7fCvGUFNmtsXMVsbp54FHCR9aZxI+QIm/q/0+GEHSHOBtwNeA5A6d\nusURW5knmtk1AGY2aGbP1jOGlFZgqqRWYCrwZK3jMLO7gWdKikerc5/eH02bNMbwCmBTan4T4Y1S\nWr45llfTx4DPS/ot8Hng4gox1cJhwMJ4mqxH0psbEAOSFgObzOzhkkV1jSPlA4TWZSNiOARIPwJd\nr2MukDSP8E3zPqDLzHrjol6gqw4hfBH4FOGLVaKecRwKbJP0dUm/lvRVSQfUOQbMbDPwBeC3hGTR\nb2bL6x1HNFqd+/T+aK1ObHtH0nJCE7vUJWb2g3rHA2PG9BngVOBjZnazpHcA1xBOkZSz1/cyV4ih\nFTjYzI6L11SWAq+qdgwZ4rgYSJ8LHev+71q8FoX/EUmfAfaY2Q21iCGDht63Lmka8D3gAjN7Thr+\nU5iZ1frZJklnAFvNbIWk7nLr1CGOVuAo4MNm9itJXwIuqnMMSDqY8A1/HvAs8B1J76l3HKUy1Jk5\nnoYmDTMb7QN3LJsJ59ITcwiZcjPDpyeS8s2M01gxSfqGmX00zn6X0BQfLaZx150xhr8Bborr/Spe\nhJ5R7RjGikPSEYRvdg/FD6g5wIOSjq12HJX+RySdRzgtckqquOqvRQWl9c2l+JtczUhqIySM68zs\nlljcK2mWmW2Jp2231jiM44EzJb0NmAxMl3RdnePYRGj5/irOf5fwxWZLnV+LU4ENZvY0gKSbgD9s\nQBww+uu/T++PiXJ6Kv0t9lbgHEntkg4lnK6538y2ADviHRMCzgVuKbOvffGkpJPi9FuB5JpL2Ziq\nXHfillg3kg4H2s1sez1jMLPVZtZlZoea2aGEN+xRsSlctzgkLSKcEllsZrtSi+r59wB4ADhM0jxJ\n7cDZMYaaiv/nVwNrzOxLqUW3Ei6+En9X+31QxMwuMbO58X/hHOBnZnZuPeOI7//fxfcEhA/vR4Af\n1CuG6DfAcZKmxL/PqcCaBsQBo7/++/b+qOYV/Gr+AH9GOE+8E9gC/Di17BLCxZu1xLuZYvnRwKq4\n7Ms1iOkEwgfESuB/gCMrxVSDGNqA6+JxPgh01zuGMjE9Qbx7qs6vxTrCm3RF/LmqUa8F4W66x2Kd\nF9fpdV9AuIawMvUaLAI6gZ8SvtQsAzrq+L9wEsN3T9U1DuCNwK+Ahwit8YMa8VoASwg3JawiXIBu\nq3UcwI2Eayh74ufm+8eqc1/eH96NiHPOucwmyukp55xzTcCThnPOucw8aTjnnMvMk4ZzzrnMPGk4\n55zLzJOGc865zDxpOOecy8yThnvJUBh7ZUXq58JY3iPp6CrX9TFJU1LzP5Q0vYr7Pncc67dLuktS\nSzXqd24s/nCfe8mQ9JyZHVim/E7gE2b26yrWtQF4s8U+hqq431bCk/5H2nC371m2+3tC9+xjddro\n3D7zlobbr0g6TdIvJT0oaamkAxRG3VuaWqdbw6PP/YfCKI2rJS2JZR8ldC99p6Q7YtlGxdELJX1c\nYdS2VZIuiGXzFEaU+0rc1+2SJpcJ8a3Ar5OEEVtJV8YY1kh6s8IohY9L+j+p7W4B3l31F8y5Ep40\n3EvJlJLTU+9IL4y9AX8GOMXMjiZ8o/84YSTEY1Onm84m9OUDoQv2txD6NTpJ0hFm9mVCPz/dZpb0\nrmuxjqOB8wiD2hwHfFDSm+I6rwH+zcyOAPqBPy9zDEn/ZgkDdscY/hP4PmEgriOA82JX3BA653sL\nztVYQ7tGd67KdprZkaMsE+FDfD7wy9ilezvwSzMbkvQTQvfe3yN0tf7JuN3Zkj5IeK/MjtuvHqOO\nBcBNZrYTCl1jn0joWXSDDQ9a9SBhzIVSswi9oqYlPeauBh6xOLCOpCeAVwLPxGPYI+kAM3thlPic\n22eeNNz+ZrmZvatM+beADwN9wANm9kLsNvoThGsXz0r6OmG8iLEYxV35i+EBbnanyoeAKYy0s0wd\nyXb5kn3kCeOTJyYB6S7inas6Pz3l9hcG3AucIOnVAPF6xmFx+V2Ekd8+yPCpqenAC4RxWroI3Z8n\nnovLS+u4G3h7HE/hAMK4zHcz9siGaY8STmONi6SXAdvNbGi82zo3Hp403EtJ6TWNy9ILLQxWdR5w\no6SHgF8Cr43L8sBthPEobotlDxHGqFgLfBO4J7W7rwA/SS6Ep+pYAfw3YVCbe4Gvxv3AyCE1y926\n+GNg4SjHZ6NsA3ByErdzteS33DrXZOJ1kAvNbP04tvke8OnxbOPc3vCWhnPN5yLCRfdMFMYJv8UT\nhqsHb2k455zLzFsazjnnMvOk4ZxzLjNPGs455zLzpOGccy4zTxrOOecy+/9gZDvYg1FFIwAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ce03b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid_id = cem.get_var_grid('land_surface__elevation')\n",
    "spacing = cem.get_grid_spacing(grid_id)\n",
    "shape = cem.get_grid_shape(grid_id)\n",
    "\n",
    "z0 = cem.get_value('land_surface__elevation').reshape(shape)\n",
    "plot_coast(spacing, z0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to elevations, the *CEM* model also takes sediment flux from *Rafem*. *CEM* can deal with many sediment sources and so takes a grid of sediment fluxes. However, *Rafem* provides only a single sediment flux value at its channel exit.\n",
    "\n",
    "We create an empty array of fluxes that we will use to pass river mouth fluxes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "qs = np.zeros_like(z0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the coupled model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're now ready for the time loop! The basic idea is to:\n",
    "1. Run *Rafem* for one time step\n",
    "2. Get its river sediment flux and elevations\n",
    "3. Pass these values to *CEM*\n",
    "4. Update *CEM*\n",
    "5. Pass *CEM*'s elevations back to *Rafem*.\n",
    "\n",
    "Note that we need to do a little work to exchange sediment discharge from *Rafem* to *CEM* because *Rafem* provides only a single value for sediment discharge and is a volume flux, not a mass flux.\n",
    "\n",
    "Because *CEM* deals with elevations on land differently than *Rafem* we also have to do a little trickery. In *CEM* all land cells have an elevation of 1., so we don't want to reset *Rafem*'s land values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set some constants for the time loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RIVER_WIDTH = dict(raf.parameters)['channel_width'] # Convert unit-width flux to flux\n",
    "RHO_SED = 2650. # Used to convert volume flux to mass flux\n",
    "N_DAYS = 150 * 365\n",
    "N_DAYS = 50000\n",
    "TIME_STEP = int(raf.get_time_step())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running rafem until 0\n",
      "running rafem until 500\n",
      "running rafem until 1000\n",
      "running rafem until 1500\n",
      "running rafem until 2000\n",
      "running rafem until 2500\n",
      "running rafem until 3000\n",
      "running rafem until 3500\n",
      "running rafem until 4000\n",
      "running rafem until 4500\n",
      "running rafem until 5000\n",
      "running rafem until 5500\n",
      "running rafem until 6000\n",
      "running rafem until 6500\n",
      "running rafem until 7000\n",
      "running rafem until 7500\n",
      "running rafem until 8000\n",
      "running rafem until 8500\n",
      "running rafem until 9000\n",
      "running rafem until 9500\n",
      "running rafem until 10000\n",
      "running rafem until 10500\n",
      "running rafem until 11000\n",
      "running rafem until 11500\n",
      "running rafem until 12000\n",
      "running rafem until 12500\n",
      "running rafem until 13000\n",
      "running rafem until 13500\n",
      "running rafem until 14000\n",
      "running rafem until 14500\n",
      "running rafem until 15000\n",
      "running rafem until 15500\n",
      "running rafem until 16000\n",
      "running rafem until 16500\n",
      "running rafem until 17000\n",
      "running rafem until 17500\n",
      "running rafem until 18000\n",
      "running rafem until 18500\n",
      "running rafem until 19000\n",
      "running rafem until 19500\n",
      "running rafem until 20000\n",
      "running rafem until 20500\n",
      "running rafem until 21000\n",
      "running rafem until 21500\n",
      "running rafem until 22000\n",
      "running rafem until 22500\n",
      "running rafem until 23000\n",
      "running rafem until 23500\n",
      "running rafem until 24000\n",
      "running rafem until 24500\n",
      "running rafem until 25000\n",
      "running rafem until 25500\n",
      "running rafem until 26000\n",
      "running rafem until 26500\n",
      "running rafem until 27000\n",
      "running rafem until 27500\n",
      "running rafem until 28000\n",
      "running rafem until 28500\n",
      "running rafem until 29000\n",
      "running rafem until 29500\n",
      "running rafem until 30000\n",
      "running rafem until 30500\n",
      "running rafem until 31000\n",
      "running rafem until 31500\n",
      "running rafem until 32000\n",
      "running rafem until 32500\n",
      "running rafem until 33000\n",
      "running rafem until 33500\n",
      "running rafem until 34000\n",
      "running rafem until 34500\n",
      "running rafem until 35000\n",
      "running rafem until 35500\n",
      "running rafem until 36000\n",
      "running rafem until 36500\n",
      "running rafem until 37000\n",
      "running rafem until 37500\n",
      "running rafem until 38000\n",
      "running rafem until 38500\n",
      "running rafem until 39000\n",
      "running rafem until 39500\n",
      "running rafem until 40000\n",
      "running rafem until 40500\n",
      "running rafem until 41000\n",
      "running rafem until 41500\n",
      "running rafem until 42000\n",
      "running rafem until 42500\n",
      "running rafem until 43000\n",
      "running rafem until 43500\n",
      "running rafem until 44000\n",
      "running rafem until 44500\n",
      "running rafem until 45000\n",
      "running rafem until 45500\n",
      "running rafem until 46000\n",
      "running rafem until 46500\n",
      "running rafem until 47000\n",
      "running rafem until 47500\n",
      "running rafem until 48000\n",
      "running rafem until 48500\n",
      "running rafem until 49000\n",
      "running rafem until 49500\n"
     ]
    }
   ],
   "source": [
    "for time in xrange(0, N_DAYS, TIME_STEP):\n",
    "    if time % 500 == 0:\n",
    "        print 'running rafem until {time}'.format(time=time)\n",
    "        \n",
    "    raf.update_until(time)\n",
    "    \n",
    "    sea_level = raf.get_value('sea_water_surface__elevation')\n",
    "    \n",
    "    # Get and set sediment flux at the river mouth\n",
    "    raf_qs = raf.get_value('channel_exit_water_sediment~bedload__volume_flow_rate')\n",
    "    \n",
    "    y, x = raf.get_value('channel_exit__y_coordinate'), raf.get_value('channel_exit__x_coordinate')\n",
    "    qs[int(y[0] / spacing[0]), int(x[0] / spacing[1])] = raf_qs[0] * RIVER_WIDTH * RHO_SED\n",
    "    \n",
    "    cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
    "\n",
    "    # Get and set elevations from Rafem to CEM\n",
    "    raf_z = raf.get_value('land_surface__elevation') - sea_level\n",
    "    cem.set_value('land_surface__elevation', raf_z)\n",
    "\n",
    "    # print 'running cem until {time}'.format(time=time)\n",
    "    cem.update_until(time)\n",
    "    \n",
    "    # Get and set elevations from CEM to Rafem\n",
    "    cem_z = cem.get_value('land_surface__elevation')\n",
    "    raf.set_value('land_surface__elevation', cem_z + sea_level)\n",
    "    \n",
    "    qs.fill(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the final elevations as seen by *Rafem*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAADLCAYAAACS0c8UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYXVWZr9/fOZWQQAhFEUyAoEEGIaLNYCAMQiGIURFw\nAhwQ1IvdzW1ERRDw2i3dz6UBLzh10+0EF0HQyNSAokSaIIpMIYEMBJJLIqawioQQwpCQqjrf/WOt\nfWrXqVNVu1JnqqrvfZ56zt5r/O1TZ+9vr+lbMjMcx3EcJwu5egtwHMdxRg5uNBzHcZzMuNFwHMdx\nMuNGw3Ecx8mMGw3HcRwnM240HMdxnMxU3WhIWi3pSUkLJT0Sw1okzZP0jKR7JDWn0l8kaYWk5ZKO\nr7Y+x3EcJzu1aGkY0GpmB5rZITHsQmCeme0D3BvPkTQTOBWYCcwBrpbkrSHHcZwGoVYPZJWcnwhc\nF4+vA06OxycBN5lZp5mtBlYCh+A4juM0BLVqafxO0mOSzophU82sIx53AFPj8a7AmlTeNcBuNdDo\nOI7jZKCpBnUcYWZ/lbQzME/S8nSkmZmkgXyZuJ8Tx3GcBqHqRsPM/ho/10q6jdDd1CFpmpm1S9oF\neCEmbwN2T2WfHsOKDGJgHMdxnH4ws9KhgiGjajoslLQtkDezVyRtB9wDXAIcB7xoZpdLuhBoNrML\n40D4jQTDshvwO2AvS4mUZGvb1tDV2Umh0A1AobubQqFQPO7q7KTQHeMK3RS6C8W06bjSvF2dWyh0\nF1L5kjIKdG3ZUjzuyddTf7rsQqHAL+/+LR9///t66ejq7OyVvlAo9Mmb1txfvnJldHeFNLl8jlw+\n/C5yORWPm8blyOVieF7F47v/9AQfPPIAmsblinEhb6533jJlJmWUxpfqKJZdpv5cXr3q7h2eL+a7\n+if38Q9feE+fevqvv6espvH5fuN61Z+Oz/eUn4SrKUdX/CmmPzenjr/1z7/iy9/4QL/xXan7rZgG\n2Jz8BlNpulJpesXH/Aseel+xrO6C0d3Pvdy2cVMxrrtgPHvbtbzl5DPLpq0E3YUeHel6B4ovF77u\n1zcw5QOfptusV1y5tIPV2V9dvcIGiSddXhJeMPj9zXDUx3rC0unL5UmXZX3j7aozqAaSKmI0qt3S\nmArcJimp62dmdo+kx4C5kj4PrAZOATCzZZLmAssI98zZ5m54q07y0HQcJyNSXyMxRqiq0TCzVcAB\nZcLXE1ob5fJcClxaTV2O4zjO1uGvmFVi5l571ltCJvZ+87R6S8jErANn1FtCJmYftXe9JQzKDvv2\neY9rSLbd+531lpCNt8yst4Ka4kajSrx9773qLSET+4wUo3HQHvWWkInDjm4so5H0zedDFzH5nGh2\no1FZ3Gg4jjNaaH91M9D/4LDjDBU3Go4zSnED4VQDNxoOXZ2FektwqkA+N+zZlY7TBzcajjOKmTZp\nQq+1Bm5InOHiRsNxxhDeZeUMl0zrNCRNBY4gOBR8HVgCPGZm3q8xCvDFfaOT7oIVB8LzUr8rxh1n\nKAxoNCQdQ9jrYifgcYKPqAkEV+Z7SfolcKWZbay2UMdxhkbaYCSf3e7/0xkmg7U0PgCcZWbPlUZI\nGgecABwP3FwFbU6N6OoseGtjjJDPie5uNxzO1jOg0TCz8weI6wRuq7gip+a4wXAcJytZxzR2BD4D\nzEjlMTP7YpV0OY5TIbrNimMaPhDuDJesr5i/Bt4CPAk8BiyIf4MiKS9poaQ743mLpHmSnpF0j6Tm\nVNqLJK2QtFzS8UO7FGdr8XUaox8fBHcqRVYvt9uY2Ve2so5zCa7Ot4/nFwLzzOwKSV+L58leGqcC\nM4l7aUjax2doVR/vnho7+JiGM1yyPi1ukPQFSbvElkKLpJbBMkmaThhM/zGQrCo6EbguHl9HmIkF\ncBJwk5l1mtlqYCVhMyanynhLY+zg3VPOcMna0ngDuAL4OpA8YQx46yD5vg2cD0xOhU01s4543EHY\nqAnCGpCHUunWEFocTpXxlsbYwVsaznDJajS+Sth2dV3WgiWdALxgZgsltZZLY2Y2yJ7f/ut2nGHy\n5h22pW3jJsDHNpzhk9VorAA2DbHsw4ETJX2AsCBwsqTrgQ5J08ysXdIuhAWDAG3A7qn802NYH664\n8koK3QXMChw2+1Bmz5o1RGmOMzbYbfJEIPigan91c3GhnzP6mT9/PvPnz694uVmNxuvAIkn3Ebqq\nYJApt2Z2MXAxgKSjga+a2emSrgDOAC6Pn7fHLHcAN0q6itAttTfwSLmyLzjvPLo6OykUugEodHdn\nvAynHL64b/TStnETu02eWFwdvqXbx6/GCq2trbS2thbPL7nkkoqUm9Vo3B7/kratGHrXUZL+MmCu\npM8Dq4FTAMxsmaS5hJlWXcDZZt6Wdpzh8tzLr9dbwuhjDD+ashqNJWb2WDpA0oeyVmJm9wP3x+P1\nwHH9pLsUuDRruU5l8FaG4wwRacwajqxPix9KekdyIukTwDeqI8lxnGrg4xlOJchqND4GXCdpX0ln\nAWcD762eLMdxhksyCJ6XigZjfN5blc7wyPQLMrNngU8QHBR+FHifmb1cTWFO7fDFfaOXaZMm1FuC\nM8oYbD+NxSVBLQRD87AkM7N3Vk2ZUzN8TGN0kt6EyXEqxWAD4ZkHux3HaSzcYDjVYDCjsc7MXh0o\ngaTtzeyVCmpyHMdxGpTB+iX+S9KVko6StF0SKGlPSZ+XdA8wp7oSHccZDu46xKkkg7U0jiN4qf07\n4PDo2bYLeBr4FfAZM2uvrkSn2viK8NFPshGTrwh3hstg270awTj8qjZynHrgBmN0stvkibRt3FTc\ntc9xKoE/LRxnFJNMufV1Gk6l8F+Q4+s0xgDTJk0oLvZznOHgRsNxRjnTJk0gn1Px2HGGQ2ajIend\nkj4bj3eWtMcg6SdIeljSIklLJH0zhrdImifpGUn3SGpO5blI0gpJyyUdv5XX5AwRH9MY3bS/utm3\neXUqRqanRXzgXwBcFIPGAzcMlMfMNgPHmNkBwAHAHEmHAhcC88xsH+DeeI6kmcCpwEzCNN6rJfnT\nzHGGSbql4TjDJetD+cPAScBrAGbWBmw/WCYzSxz5jwfGEfbUOBG4LoZfB5wcj08CbjKzTjNbDawE\nDsmoz3GcMrixcCpNVqPxhpkVR0vTC/0GQlJO0iKgA7jHzB4BpppZR0zSAUyNx7sCa1LZ1xB28HMc\npx+Wrd045Dw+g8oZDlk3YfqlpB8AzZK+AHwO+PFgmaKhOUDSDsBtkvYviTdJA3W2ekes46S49s7D\nUmcbmbnzZJat3cjbdurb8O9vHMMX+DnDYVCjIUnAL4B9gVeAfYBvmNm8rJWY2ctxf/H3AR2SpplZ\nu6RdgBdisjZg91S26TGsD1dceSWF7gJmBQ6bfSizZ83KKsVJ0TJtWzas3VRvGQC869jp5MfnWfC7\nv9RNwwm37s9dH1lSt/rL8f3be/fQfvZDf6IrHi946H0ARcOxV8ukPvkTp4U+rjH2mD9/PvPnz694\nuVlbGr82s/2Be7IWLGkK0GVmGyRNJGzadBlwB3AGcHn8vD1muQO4UdJVhG6pvYFHypV9wXnn0dXZ\nSaHQDUChuzurLCfSvPNE1re/Tsu0bcmlHiavvvzGoPkgzLiqlMGZechUHrt3Dbm82P2YXWi7v/ae\naU68/R3c9ZElnHBraAz/5rRlAMz5+cximnmfebqqGv7lloO56COP8a3bZnH+hx/l27fN4pyTwy2w\nuRBaB139tL3fttP2PP1i8Btazni4wRh7tLa20traWjy/5JJLKlLuoEYjdiEtkHRIHJPIyi6E3f7y\nhLGTX5jZryU9BMyV9HlgNXBKrGeZpLnAMoJ/q7OjG5Mxy0svtLPDlJ2L551bwryC5CGfywsKeSZO\nyvHqhvBA37h+4Ic+wOSWCcUH/oa1m3pNuQ1llHep3TxlIhvWbgr1AlN2CUNb5dJP2XVb1nf0NSq7\nvnUy7X/ucYq878FvYtkjHcXztvvb2eaIKXQ+9OKg11EpTrz9Hdxx8mJA3PWRJeTyYs7P3w4E45F8\n3+/96X4A3H36cgA+dMN+3PapZRXRcOmt7+IbH32Mf7nlXZz/4Uf51m2z+PKHH+3XSJSjnLFI6C5Y\n0XCMz+e8i8rZapTluSzpaWAv4M/EGVQEe1LzTZgk2dq2NX1aGoX4Jlbo7g5xsfVRKHRT6C4U06bj\nSvN2dW6hEG+mkC8po0DXli3F4558vVs6SdlJmnRcoRB09Q4r9MmbpN2wdi2TW3bi5RdD710uL8aN\n37ZPueO26WbTqwW2m5xn4/rNNO88kVdfDlpzORUf8E3jcuRyKhqFtOFJG41cXjRPCcYnnbd5ygQ2\nvvRGrzKTMqbsuh0b12+mZeq2xfTrOzbxpunbse6vrzN9zx2KZbf/+RVm7Lcja1a+zF7v3Ilnl64v\n1p/Lq1jm2lmT2W1xXyOZy4tznzqK77/9gaLmpnG5PprSZTWNz/cb97HfHBANRt96it9JDFdTji4z\n3n/9vgDc+emn+MBP92Xup4Ph6DKjK95P6c/NqeOu1P2WHF9880Fc/JHH+paRStPT0rA+3VMQjEKp\nf6n0JkzJavDugrFqw2tUk/RYSqKpV1iZ+AHDzfqMz5SmHazO/uoqrWegeErHiMx6h5U+SwsDxCfh\n1jferjqDaiAJMxt2kzNr91Ty60yu0Nu6VWTj+vVMbtkJgO133KmPsUnzxibFPJvj5xvFB3yayS3b\nkMvlerUUEkq93G5Yt6nYitiwbjOTW7bpU16aDWs30TJ1W9Z3vJ568OZY99fXedP07Wh/LrQskrjn\nntnAW9/ewsonX+yjJWHq469w9WsncXbzXZz3/1oB+PY+93PuU0fx3f1+zzlLexuOreFjvzmAWz/4\nBEN1jHD36cuLD/M7Tn+KE68PXVg3fmrpgPm+MveA4vEVH18IBIPxzx97nM1VevGfNmmCb8ZUDcZw\nJ0jWPcJXA82ENRYfAnaIYc4ooJorwl9YU5232u+//QHOWfruqpRdLa46ZVHRWDgjHI3d9+asK8LP\nJawA35mwruIGSV+sprCxyoZ1a5nc0jKsMja+FFobCS1TQ3dTf2MV5VjfsYkN6zYX8w7GhnVDGxRf\n/dRLg6Y5u/kuzl91DFfuOZ8r95zPl585mu/u9/ti/NYYjtPuPQgIrYyb5ywaUt7+mPvpZcz99DI+\n+bO395vmi7/4m+LxFR9fyAW/PJCLbz6ISz/2eEU0lCOfE/mc3N+UU1GyjmksBmab2WvxfDvgITN7\nR5X1ldMyYsc01ne09yp/cstOvfJuWPtCCIuak3TluqdKx0W6u0KaXL6nf795yoQwy2nd5uKYRkij\nAcc0AHK5XO/xkNSYQen4Qfkxhd75e8eVr79veJmxiJJ6AM5Z+m5+cPCDxfCzHg1rGa49/GEAPvPH\nnmmrPz/2cU6796CiwehVf8YxDeh/zOKTP3s7P/3kkl7xX/j5O/neqU/0O6bRRcmYRQXHNPI59RoE\nLx3TGJ/PkZfY1FW5GYg+ppFKP4bHNAAK/Rw7JWxc/2JxTAJgfUeYQpqEJQ/8DevWMmmH5rJ5KqPj\njX7HDNKM9J37/uOAP/K3C47gBwc/yFmPHsaPZv2JXF589sHZdJtx4zELgJ6H/9zjF1KtYbmffnIJ\nn7lxf675RBhcTwxGPehvmm2yKVNiMAAmNuXpNivOqprYlC+mr6RBcUY+WY3GtcDDkm4l3G0nA9dU\nTdUIZ3LLTmxc/2KxhdAydVqvFkM6XWh99G511JqRbDASfnDwg/ztgsP50aw/FcOSlgbky2eqEsFw\nhEb4D097ks11GjNt27ipbNdUf25E8lIvYwE9A+luOJyETEbDzK6SdD9wJGEG1ZlmVtcRvY3rXyx2\n8TQSG19az6Qdmvu0KvqjOa7DKDUotWSktzQSQhdVYwxQJi2NoayzqDTJKvDugvXqohqMZJvYZJru\ntEkTqj5N1xk5DKV7ahHQHvOYpDeb2XPVkTUw6zvai33/YT1D34HjYFQKfeI2rn+x+GZfaTasW8vk\nHVt6jWmMBEaDwXD6kjYS7a9u7rVz326TJ/Zax5EOB9+syemfTEZD0jnAPxH8RKVfnWs+EA6huyd5\nM2/eeefimEFXZyfNU3buNT6w7q9tNE95ExvWru01+Lxh3dpiecm4AvS0YJLw9R3tNO/8pt7xZYzO\n+o52mqfsPCJdmoyWlobTP8W9wnMqGoZkZlXacAylReKMTbK2NL4EvM3MaufbYQg0T9m5+Ha/7vm2\nXuMDLVOnxQd/726gpFuo0N1dNAzJKmxIVmW/QMvUaax7vi2GFYrx6zv+2qu8tCEbabjBGP2kDcFA\nhsENhjMYWY3Gc8DQHffXgXIP7+Ypbyq2HsrRvPObioalkPLJkxiIlmm79Jlymy5zJLYu0nhLY3SS\nGIf0dFsIXVX9eb31loYzGAMaDUnnxcNngfmS7gK2xDAzs6uqKa6WBGOzZfCEjjOCKDUCSVdUYjhK\nxzSScMfpj8FeL7cHJhFaGvMI27ZOiuGDbvcqaXdJ90laKmlJsopcUoukeZKekXSPpOZUnoskrZC0\nXNLxW3thTna8lTF26DYr/rW9sqnPgkAIhqP91c3FWVc+c8pJM2BLw8y+WRoWXZ1PMrOXM5TfCXzZ\nzBZJmgQskDQP+Cwwz8yukPQ14ELgQkkzgVOBmYQ9NX4naZ/0VrOO42Sjv26m9Orx8fkcm7q6e63d\nSBb8+foMpxxZfU/dKGlydB+yGFgm6YLB8plZu5ktisevAk8RjMGJwHUx2XWExYIAJwE3mVlndIi4\nEui9dZlTcbo63SaPRpKWQvKXtB7SJCvA05+JoXCD4ZQja7/E281sI+HhfjcwAzh9KBVJmgEcCDwM\nTDWzZOedDoITRIBdgTWpbGsIRsapIt49NXpJupraX91cNAzluqSg997hbjCc/sj6tGiSNI5gNO40\ns0569tYYlNg1dQtwrpm9ko6Lu/MNVNbYdVzvOMMgPcg9bdKEft2HOM5QyDrl9geErVmfBH4fWw1Z\nxjSIxuYW4HozS/YD75A0zczaJe1CWDQI0Absnso+PYb14oorr6TQXcCswGGzD2X2rFkZL8Nxxh7J\nYr7dJk/kuZdfp7vb38PGAvPnz2f+/PkVLzer76nvAd9LziX9GThmsHySBPwEWGZm30lF3QGcAVwe\nP29Phd8o6SpCt9TeQJ99yS8477w+rtEdxxmY0vEMZ3TT2tpKa2tr8fySSy6pSLlD8T1VJHYpdQ2a\nEI4APg08KSlxcHgRcBkwV9LnCS2YU2K5yyTNBZbF8s+2LBt+OMPCF/eNbhLng+UGwh1nqGyV0ciK\nmf2B/sdNjusnz6XApVUT5fTBDcbop23jpuL6DMcZDv60cJwxgBsLp1JkXadxiqTJ8fgbkm6TdFB1\npTm1wtdpjG7S258mC/ecYTKGjXDWlsY3zGyjpCOBYwmD2/9RPVlOLfHuKccZImPY+GZ9WiTTk04A\nfmRmdxH8UDmO08CkN15ynEqQ1Wi0SfohwS/UryRNGEJex3HqQHF9xvY9my45znDJ+uA/BfgtcLyZ\nbQB2BM6vmirHcYaNT691qkFWozEN+JWZrZB0DMGI9Fl054xMfCB8dJKsy2h/dXNxANxbG85wyWo0\nbgW6JO1FcCkyHbixaqqcmuID4aOXxP+UT7l1KkXWp0XBzLqAjwDfN7PzgV2qJ8upJd7SGJ2kB8GT\nloZ3WTnDJavR2CLpk8BngLti2LjqSHJqjbc0Ri/prVvdYDiVIOvT4nPAYcD/NrNVkt4KXF89WY7j\nVJJuMx/PcCpCJqNhZkuBrwJLJO0P/MXMLq+qMqdmePeU4zhZyepGpBV4Bvh34GpghaSjM+S7RlKH\npMWpsBZJ8yQ9I+keSc2puIskrZC0XNLxQ74ax3F6kQyE++wpp1Jk7Z66irBG4ygzOwo4Hvh2hnzX\nAnNKwi4E5pnZPsC98RxJMwmLB2fGPFdL8s72GuBjGqOTto2bisdp/1OOMxwyb/dqZk8nJ2b2DBnc\nqpvZA8BLJcEnAtfF4+sIW8gCnATcZGadZrYaWAkcklGf4zj94NNtnUqSdT+NBZJ+DNwACPgU8NhW\n1jnVzDricQcwNR7vCjyUSreGsHuf4ziO0yBkbWn8HfAU8EXgHGAp8PfDrTzuyjfQa5C/IjnOVpKs\n08hLPX8+puEMk0FbGpKagCfMbF/gygrU2SFpmpm1S9oFeCGGtwG7p9JNj2F9uOLKKyl0FzArcNjs\nQ5k9a1YFZDnO6GPapAnFwXDA99MYQ8yfP5/58+dXvNws4xJdkp6W9BYz+3MF6rwDOAO4PH7engq/\nUdJVhG6pvenHv9UF551HV2cnhULw2F7o7i6XzHHGPGmDAT6+MZZobW2ltbW1eH7JJZdUpNysYxot\nwFJJjwCvxTAzsxMHyiTpJuBoYIqkvwD/CFwGzJX0eWA1wfkhZrZM0lxgGdAFnB27rxzH2Qp89pRT\nDbIajW+UnGf65ZnZJ/qJOq6f9JcCl2bU5DhOBrx14VSSAY2GpL0Js53ml4QfCfy1irqcGtLVWfC1\nGqOY3bafGNyku/FwKsBgT4rvABvLhG+Mcc4owA3G6GW3yRPJ58S0SRN89pRTEQZ7Wkw1sydLA2PY\nHtWR5NQa9z01emnbuKm4EZPjVILBjEbzAHETBohzRhDe0hjdtL0SBsS9e6qCjOHvcrCnxWOSvlAa\nKOksYEF1JDmO4zQ4Y3i9y2Czp74E3CbpU/QYiYOBbYAPV1OY4ziVIS95K8OpGAMajbhq+3DgGGB/\nwlTbu8zsv2shznGcyuLrNJzhkmVFuAH/Hf8cxxkh7DZ5Yq8Ffo5TCXwE1PHZU6OY9B7hjlMJ3Gg4\nPntqDOA79zmVwp8WjjNKKbc+w8c0nOHiRsNxRilpg+Gzp5xK0XBGQ9IcScslrZD0tXrrcRzHcXpo\nKKMhKQ/8GzAHmAl8QtJ+9VW1dSxdsbLeEjLxzHPt9ZaQiUcfX1VvCZn40/0r6i1hQPI5sWH5onrL\nyMTrK/p4MGpM/rys3gpqSkMZDeAQYKWZrTazTuDnwEl11rRVLFv5/+otIRMrRorRWLi63hIy8dDv\nG9NoJNu9ArzsRqOyuNGoK7sBf0mdr4lhjuNUAN/u1RkujWY0fLTOcSpEuTUaPiDuDBc10o6qkmYD\n3zSzOfH8IqBgZpen0jSOYMdxnBGEmQ27qdloRqMJeBo4FngeeAT4hJk9VVdhjuM4DpB9j/CaYGZd\nkv4B+C2QB37iBsNxHKdxaKiWhuM4jtPYNNpA+IA0ysI/SbtLuk/SUklLJH0xhrdImifpGUn3SGpO\n5bko6l4u6fga681LWijpzkbVKalZ0s2SnpK0TNKhDarzy/F/vljSjZK2aQSdkq6R1CFpcSpsyLok\nHRyvbYWk79ZA47fi//wJSbdK2qGeGvvTmYo7T1JBUkuj6pR0TvxOl0hKjwdXRqeZjYg/QnfVSmAG\nMA5YBOxXJy3TgAPi8STCOMx+wBXABTH8a8Bl8Xhm1Dsu6l8J5Gqo9yvAz4A74nnD6QSuAz4Xj5uA\nHRpNJ2H697PANvH8F8AZjaATeDdwILA4FTYUXUmvwyPAIfH418CcKmt8b/KdAJfVW2N/OmP47sBv\ngFVASyPqJOx9NA8YF893rrTOkdTSaJiFf2bWbmaL4vGrwFOEB8qJhIcf8fPkeHwScJOZdZrZasI/\n7JBaaJU0HfgA8GMgmTnRUDrj2+W7zewaCGNbZvZyo+mMNAHbxkkb2xImbNRdp5k9ALxUEjwUXYdK\n2gXY3sweiel+mspTFY1mNs/MEt/8DwPT66mxP52Rq4ALSsIaTeffA/8an5GY2dpK6xxJRqMhF/5J\nmkGw9g8DU82sI0Z1AFPj8a4EvQm11P5t4HwgvWlGo+ncA1gr6VpJj0v6kaTtGk2nmbUBVwLPEYzF\nBjOb12g6UwxVV2l4G7XV+znCmy5ltNRVo6STgDVmVrpMvaF0AnsDR0l6SNJ8Se+qtM6RZDQabsRe\n0iTgFuBcM3slHWehrTeQ5qpfj6QTgBfMbCE9rYzeIhpAJ+Ht/SDgajM7CHgNuLCXiAbQKWlHwtv7\nDMLNNknSp3uJaACdZSsdXFddkfR1YIuZ3VhvLaVI2ha4GPindHCd5AxGE7Cjmc0mvCzOrXQFI8lo\ntBH6FBN2p7eFrCmSxhEMxvVmdnsM7pA0LcbvArwQw0u1T49h1eZw4ERJq4CbgPdIur4Bda4hvMU9\nGs9vJhiR9gbTeRywysxeNLMu4FbgsAbUmTCU//OaGD69JLzqeiWdSehC/VQquJE07kl4UXgi3kvT\ngQWSpjaYTmLdtwLE+6kgaUoldY4ko/EYsLekGZLGA6cCd9RDiCQBPwGWmdl3UlF3EAZGiZ+3p8JP\nkzRe0h6EJuQjVBkzu9jMdjezPYDTgP82s9MbUGc78BdJ+8Sg44ClwJ2NpBP4MzBb0sT4GzgOWNaA\nOhOG9H+O/4eNCjPXBJyeylMVJM0hvBGfZGbpHaMaRqOZLTazqWa2R7yX1gAHxa6/htEZuR14D0C8\nn8ab2bqK6qzkaH61/4D3E2YqrQQuqqOOIwljBIuAhfFvDtAC/A54BrgHaE7luTjqXg68rw6aj6Zn\n9lTD6QT+BngUeILwprRDg+r8JmHiw2LC4PK4RtBJaEk+D2whjP19dmt0AQfHa1sJfK/KGj8HrCAY\n4+Q+urqeGkt0vpF8lyXxzxJnTzWazvh7vD7WuwBorbROX9znOI7jZGYkdU85juM4dcaNhuM4jpMZ\nNxqO4zhOZtxoOI7jOJlxo+E4juNkxo2G4ziOkxk3Gs6IQdLJ0S3121JhM8q5sK6BllZFV/M1qOs7\nko6Mx6vTbrmHWM4HJV1SWXXOWMONhjOS+ATwQPwc0UjKdO9J2gk41Mz+EIO2emGVmf0K+JCkiVtb\nhuO40XBGBNE55BHA/yC4RCmXZkL0lPtk9JbbGsPPVNjg526FDYnSG9N8XtLTkh6O3nW/X6bcoxU2\nsVoYy50UoyZJ+mXc8OaGVPpjY7onJf0kur1JWgmXSVoAfFzS8ZIelLRA0tzo2beUjwJ3l9E0MV7P\n5yW9RWE8A5ihAAANdklEQVRjnWvjtdwg6ThJf4jXOyuVdT5wwoBftuMMgBsNZ6RwEnC3ma0AXpR0\nUJk0/xPoNrN3Eloj10naJsb9DXAK8A7gVEm7SdoV+F/AoQSD9DbKv8mfB5xtZgcSXMhsiuEHAucS\nNrh5q6TDJU0ArgVOiTqaCHscEMteZ2YHA/cCXweOjecLCJtllXJ4jEuzPcGX0M/M7CcEj6t7Av8H\n2Df+nWZmRwJfJbiPSHiMsHmP42wVbjSckcInCDvlET/LdVEdAdwAYGZPE3wa7UN4WN9rZq+Y2RsE\nJ4MzCBsi3W9mGyx4rf0l5V1e/xH4tqRzCG6nu2P4I2b2vAVfPIsI+4K8jeAJd2VMcx1wVKqs5Bpm\nE4zNg5IWAp8B3lym7l2AtalzAf8FXGNmN6TCV5nZ0qhlKcEoASyJ15qwluDW3XG2iqZ6C3CcwYgD\nv8cA+0sywta/RvCO2id5P8W8kTruJvz2S1sV/e05crmku4APAn+U9L4hlpkOey11PM/MPtmP3oRN\nQHoMwoA/EJx33pQKT2spEJwCJsfp+3wCPS0lxxky3tJwRgIfA35qZjMsuKd+M7BKUmk3ywPEPRmi\nW+g3Ezx6ljMGRvCqe7SkZoXtWz9Kme4pSXvGt/grYp7+urGM4IV5hqQ9Y9jpwP1l0j4MHJGkk7Sd\npL3LpHsK2Ksk7B+BlyT9e5n0g7EPwaOp42wVbjSckcBpwG0lYbfE8PSOdFcDOUlPEvaQP8PCXsll\nd60zs+eBSwl7XPwBWAVsLFP/uZIWS3qC8AafDEyXK/MNgovqX0YdXcB/lqa3sHfzmcBNsdwHCcao\nlF8BrekqYv5zgYmSLuvn+qyf49ZYpuNsFe4a3RnTSNrOzF6LLY1bgZ+Y2X/VW1caSQ8AJ5jZy8Ms\nZyph8Py4yihzxiJuNJwxjaRvEXbgmwD81sy+VGdJfZB0CLDJzIbVrSTpXYR9uJ+sjDJnLOJGw3Ec\nx8mMj2k4juM4mXGj4TiO42TGjYbjOI6TmRG3uC8u7nIcx3GGiJn1t/g1MyPOaABwzDdBqWuXIFdy\nDj1hyWc6T059z8vlzZInSVeat5y2/vKUll+aZ6B6y117f+FJ/oG+rzIa8hL5kvNen7ne50nYYOfl\n8m5VnkpoLZNnQA0l9Rw8+7fFG6pJYkIuR1PM1xTDkvMJuVyfsCaJCSXnvfKUnDdJWFeBQsEodId3\nqUK3USgYXZ2FXufF+ILRtaWbQsH6xBfzxPNeeTq7e533Sp8qIym3q7MQyyiUz1OSPokLYYU+YeXy\n9KRPXUt3gXzTOHL5PAC5XI5cPl88bxo3jlwuTy6fi/EhrmncuHCezxfzFMNimlwuV0zTNC5dR+/4\npvHje+rNxTT5HE3jxvcqo28duVT6njp65cmVaC2nK1VGogNA6WfAMPDuKcdxHCczbjQcx3GczLjR\ncBzHcTLjRsNxHMfJjBsNx3EcJzNuNBzHcZzMuNFwHMdxMuNGw3Ecx8mMGw3HcRwnM240HMdxnMy4\n0XAcx3Ey40Zja3lhZb0VwJrl9VYAQNezS+otgQ3LF9VbAgDL//hsvSXw+wdW1FsCTz71l3pLAODZ\n59vrLYGFi+t/f1QSNxpbixuNIt2r6n9TvNwgRuPpP66qtwQe+IMbjYRn/+pGo9K40XAcx3Ey40bD\ncRzHyYzMRtaeRr4Jk+M4ztZRiU2YRpzRcBzHceqHd085juM4mXGj4TiO42SmYY2GpI9LWiqpW9JB\nJXEXSVohabmk41PhB0taHOO+WwVNB0h6SNJCSY9KmjWYpmog6RxJT0laIunyemhI1XmepIKkllrr\nkPSt+D08IelWSTvUWkOqvjmxrhWSvlbt+mKdu0u6L94nSyR9MYa3SJon6RlJ90hqrpGefLw37qyH\nDknNkm6Ov4llkg6tx3ch6cvx/7FY0o2Stqm2DknXSOqQtDgV1m+dw7o/zKwh/4B9gX2A+4CDUuEz\ngUXAOGAGsJKesZlHgEPi8a+BORXWdA/wvnj8fuC+ATTlqvS9HAPMA8bF851rrSGlZXfgN8AqoKUO\n38V7k7KBy4DL6vFdAPlYx4xY5yJgv2p+97HeacAB8XgS8DSwH3AFcEEM/1ryvdRAz1eAnwF3xPOa\n6gCuAz4Xj5uAHeqgYTfgWWCbeP4L4Ixq6wDeDRwILE6Fla1zuPdHw7Y0zGy5mT1TJuok4CYz6zSz\n1YQLPlTSLsD2ZvZITPdT4OQKyyoQfogAzUDbAJoOqXDdCX8P/KuZdQKY2do6aEi4CrigJKxmOsxs\nnpkV4unDwPRaa4gcAqw0s9Xx//LzqKGqmFm7mS2Kx68CTxEeWicSHqDEz0rfB32QNB34APBjIJmh\nUzMdsZX5bjO7BsDMuszs5VpqSNEEbCupCdgWeL7aOszsAeClkuD+6hzW/dGwRmMAdgXWpM7XEG6U\n0vC2GF5JvgR8S9JzwLeAiwbRVA32Bo6K3WTzJb2rDhqQdBKwxsyeLImqqY4UnyO0LuuhYTcgvQS6\nVtdcRNIMwpvmw8BUM+uIUR3A1BpI+DZwPuHFKqGWOvYA1kq6VtLjkn4kabsaa8DM2oArgecIxmKD\nmc2rtY5If3UO6/5oqoy2rUPSPEITu5SLzezOWuuBATV9HTgO+JKZ3Sbp48A1hC6Scmz1XOZBNDQB\nO5rZ7DimMhd4a6U1ZNBxEZDuCx1o/nc1vovib0TS14EtZnZjNTRkoK7z1iVNAm4BzjWzV6Sef4WZ\nWbXXNkk6AXjBzBZKai2XpgY6moCDgH8ws0clfQe4sMYakLQj4Q1/BvAy8EtJn661jlIy1JlZT12N\nhpn198AdiDZCX3rCdIKlbKOneyIJb2OIDKRJ0k/N7Ivx9GZCU7w/TUOuO6OGvwdujekejYPQUyqt\nYSAdkvYnvNk9ER9Q04EFkg6ttI7BfiOSziR0ixybCq74dzEIpfXtTu83uaohaRzBYFxvZrfH4A5J\n08ysPXbbvlBlGYcDJ0r6ADABmCzp+hrrWENo+T4az28mvNi01/i7OA5YZWYvAki6FTisDjqg/+9/\nWPfHSOmeSr/F3gGcJmm8pD0I3TWPmFk7sDHOmBBwOnB7mbKGw/OSjo7H7wGSMZeymipcd8LtsW4k\n7QOMN7N1tdRgZkvMbKqZ7WFmexBu2INiU7hmOiTNIXSJnGRmm1NRtfx/ADwG7C1phqTxwKlRQ1WJ\nv/OfAMvM7DupqDsIg6/Ez0rfB70ws4vNbPf4WzgN+G8zO72WOuL9/5d4T0B4eC8F7qyVhsifgdmS\nJsb/z3HAsjrogP6//+HdH5Ucwa/kH/BhQj/xJqAduDsVdzFh8GY5cTZTDD8YWBzjvlcFTUcQHhCL\ngD8BBw6mqQoaxgHXx+tcALTWWkMZTc8SZ0/V+LtYQbhJF8a/q+v1XRBm0z0d67yoRt/7kYQxhEWp\n72AO0AL8jvBScw/QXMPfwtH0zJ6qqQ7gb4BHgScIrfEd6vFdAN8kTEpYTBiAHldtHcBNhDGULfG5\n+dmB6hzO/eFuRBzHcZzMjJTuKcdxHKcBcKPhOI7jZMaNhuM4jpMZNxqO4zhOZtxoOI7jOJlxo+E4\njuNkxo2G4ziOkxk3Gs6oQWHvlYWpvwti+HxJB1e4ri9Jmpg6/5WkyRUs+/QhpB8v6X5J+UrU7zgD\n4Yv7nFGDpFfMbPsy4fcB55nZ4xWsaxXwLos+hipYbhNhpf+B1uP2PUu+fyS4Zx/IaaPjDBtvaThj\nCknHS3pQ0gJJcyVtp7Dr3txUmlb17D73Hwq7NC6R9M0Y9kWCe+n7JN0bw1Yr7l4o6SsKu7YtlnRu\nDJuhsKPcD2NZv5U0oYzE9wCPJwYjtpKuihqWSXqXwi6Fz0j6l1S+24FPVfwLc5wS3Gg4o4mJJd1T\nH09HRm/AXweONbODCW/0XyHshHhoqrvpVIIvHwgu2GcR/BodLWl/M/sewc9Pq5kl3nUt1nEwcCZh\nU5vZwFmSDohp9gL+zcz2BzYAHy1zDYl/swQD3oga/hP4L8JGXPsDZ0ZX3BCc883CcapMXV2jO06F\n2WRmB/YTJ8JDfCbwYHTpPh540My6Jf2G4N77FoKr9a/GfKdKOotwr+wS8y8ZoI4jgVvNbBMUXWO/\nm+BZdJX1bFq1gLDnQinTCF5R0yQec5cASy1urCPpWeDNwEvxGrZI2s7MXutHn+MMGzcazlhjnpl9\nskz4z4F/ANYDj5nZa9Ft9HmEsYuXJV1L2C9iIIzervxFzwY3b6TCu4GJ9GVTmTqSfIWSMgqE/ckT\ntgHSLuIdp+J495QzVjDgIeAISXsCxPGMvWP8/YSd386ip2tqMvAaYZ+WqQT35wmvxPjSOh4ATo77\nKWxH2Jf5AQbe2TDNU4RurCEhaSdgnZl1DzWv4wwFNxrOaKJ0TOPSdKSFzarOBG6S9ATwIPC2GFcA\n7iLsR3FXDHuCsEfFcuBnwB9Sxf0Q+E0yEJ6qYyHwfwmb2jwE/CiWA3231Cw3dfFu4Kh+rs/6yQNw\nTKLbcaqJT7l1nAYjjoNcYGYrh5DnFuBrQ8njOFuDtzQcp/G4kDDongmFfcJvd4Ph1AJvaTiO4ziZ\n8ZaG4ziOkxk3Go7jOE5m3Gg4juM4mXGj4TiO42TGjYbjOI6Tmf8PuXnDHw06e6UAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d869c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = cem.get_value('land_surface__elevation').reshape(shape)\n",
    "plot_coast(spacing, z - sea_level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "z = raf.get_value('channel_centerline__elevation')\n",
    "sea_level = raf.get_value('sea_water_surface__elevation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10e16a810>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG5pJREFUeJzt3XuQVOWd//H3VxABuRoVRi6Cq0RABHFVYhFpboKoCGrw\nCgN42Y2uq6noBrO1P7EqZWU3cWNW466rgCMR0KhcQhAYDO0trKigKIiAERQMEEEQ1ARwvr8/nh5t\nJoDDdPec7nM+r6op+3b6PE9Jfc7Tz3ku5u6IiEgyHBF1AUREpP4o9EVEEkShLyKSIAp9EZEEUeiL\niCSIQl9EJEFqFfpmNtnMtpjZW1mvHWNmlWa2xswWmlmrrPfuNLO1ZrbazM4vRMFFROTw1balPwUY\nWuO1CUClu3cBnss8x8y6AVcA3TLHPGhm+kUhIlIEahXG7v4i8EmNl4cDFZnHFcCIzONLgOnuvtfd\n1wPrgLNzL6qIiOQqlxZ4G3ffknm8BWiTeXwCsDHrcxuBdjmcR0RE8iQv3S4e1nI41HoOWutBRKQI\nNMzh2C1m1tbdN5tZGbA18/omoEPW59pnXtuPmelCICJSB+5udT02l5b+HKA887gcmJX1+pVm1sjM\nOgOnAEsP9AXuHtu/u+66K/IyqH6qXxLrF+e6uefeVq5VS9/MpgP9gGPN7EPg/wE/BZ40s+uA9cCo\nTJCvMrMngVXAPuAmz0dJRUQkZ7UKfXe/6iBvDTrI5+8B7qlroUREpDA0fr5AUqlU1EUoKNWvtMW5\nfnGuWz5YVD0vZqZeHxGRw2RmeEQ3ckVEpMQo9EVEEkShLyKSIAp9EZEEUeiLiCSIQl9EJEEU+iIi\nCRJp6H/wQZRnFxFJnkhD/4wz4IorYMmSKEshIpIckYb+++/DuefCNddAnz4wYwbs3RtliURE4q0o\nlmH48kuYMwfuuw9WroRhw+CSS8J/mzSJpHgiIkUp12UYiiL0s334YbgAzJwJK1bADTfA978P7dtH\nUEgRkSITu7V3OnSAm2+GRYvgpZdg1y44/XS48kr4wx9Aa7SJiNRd0bX0D+TTT2HKFLj/fmjeHC67\nDC6+OFwMrM7XOxGR0hO77p1DqaqCdDp0//z2t+Gm70UXhQvAgAFw1FGFKauISLFIVOhnc4d33oG5\nc8NFYO1auO46+Md/hI4d81hQEZEiktjQr2nNGnjwQZg6Fc47D266CQYNUvePiMSLQr+G3bvh8cfh\ngQdC98/3vw/l5dCqVd5PJSJS7xT6B+EOL74I//3fMH8+XHppGPs/cCAcfXTBTisiUlAK/VrYvBmm\nTYPf/Q6WLoX+/UPf/5Ah0KBBvRRBRCQvFPqH6dNP4amnwi+Ajz+GG2+E8eOhTZt6L4qIyGGL3eSs\nQmvRIoT8q6/Cb34D770Hp54aJn+l05r8JSLxlriW/oHs2AGPPQYPPRTWAbrxRhgxAk46KeqSiYjs\nT907eeQeln6YPBnmzQsjfoYPD+P/Tz016tKJiCj0C6aqCpYvD/3/U6bAt78dhn5ecAGUlUVdOhFJ\nKoV+PdizJ8z6ffLJsBBcx44wahSMG6cLgIjUr8hD38x+AFwHOPAWMA44GngCOBFYD4xy9x01jiuZ\n0M+2b1/Y6Wvq1HAjOJUKC8ANGQLHHRd16UQk7iINfTNrB7wIdHX3v5rZE8A8oDvwsbv/h5n9CGjt\n7hNqHFuSoZ9t167Q+p87FxYvhi5dYOzYsBNYy5ZRl05E4qgYhmw2BJqaWUOgKfARMByoyLxfAYzI\nw3mKTvPm4SbvzJmwdSv85Cdh2OeJJ4bwf/llDQEVkeKSj+6dW4GfAF8AC9x9tJl94u6tM+8bsL36\nedZxJd/SP5itW8MQ0IcfhoYN4frrYfRoOPbYqEsmIqUu0pa+mbUmtOo7AScAzczs2uzPZJI9nul+\nEMcfD7ffDqtXh5U/ly2Dk08OE8Ceey6MDBIRiULDHI8fBLzv7tsAzOwZ4DvAZjNr6+6bzawM2Hqg\ngydOnPjV41QqRSqVyrE4xcUM+vULf598Elb//OEPw1IQ48eHCWDdu2v5ZxE5uHQ6TTqdztv35Xoj\n92xgMnAW8BfgUWApYdTONnf/dzObALSK443cunCH118PY//nzQvLP59/fhj9M3CguoBE5NCKYcjm\nROAKYB+wDLgeaA48CXQkZkM288kd1q2DBQvC3/PPQ48e4VfAqFHhRrGISLbIQ7/OJ1bo/409e0L4\nT5oULgCXXhpuAvfpoy4gEQkU+jH1pz+FEUCTJoURQNddF0YAHX981CUTkSgp9GOuegewSZNg9mwY\nMCAsAjd4MLRrF3XpRKS+KfQTpHoDmPnzw9DPsrJwE/j888Nm8E2bRl1CESk0hX5CffllGAVUWQkL\nF8KKFXD55aEb6JxzdA9AJK4U+gLARx+FewCTJ8Nf/xpa/4MHw7Bh0KxZ1KUTkXxR6Mt+3MNM4MpK\nePZZ+L//C6OAxo2D73xHG8GLlDqFvhzS5s3hF8DUqeHXwIAB4RfA4MHQuXPUpRORw6XQl1rbtCls\nArNwYfhvixZfXwD69w/bQ4pIcVPoS51UVYWbv5WV4W/JktD9U70mUOPGUZdQRA5EoS958cUXMGtW\nuBG8fDlccUW4APTurZFAIsVEoS95t2EDVFTAo4+GkT/jx4fdwLQdpEj0FPpSMFVVYQ2gKVPCxvAD\nB4ZRQEOHhqUhRKT+KfSlXuzcGfYDnjIF1q8PE8EGDw4bw2s1UJH6o9CXerd6dej/r6yEpUuhVy8Y\nNChcBPr0gSPysfOyiByQQl8i9fnn8NJLX08G270bysthzBj4u7+LunQi8aPQl6KyfHnoApoxA1q3\nDr8ALrwwLAuh+wAiuVPoS1HKngcwcya8/z5ce23YEax3by0HIVJXCn0pCe++G4aAzp4NW7aE5SCq\n7wOcdFLUpRMpHQp9KTnVy0FU/zVpEi4AgwaFYaHf+lbUJRQpXgp9KWnusHJl2BRm0aKwS9h558HY\nsXDRRdCoUdQlFCkuCn2Jld27w+5gU6bAqlVw9dXhAnDGGVGXTKQ4KPQltt57LywLXVERRv5UdwEN\nHaqNYSS5FPoSe1VV8PbboftnwYIwIWzEiPAL4Lvf1WQwSRaFviTO5s3w+OOhC+jzz8NksPJy6NQp\n6pKJFJ5CXxLLHZYtC0NBp0+HHj1C6/+yy9T9I/Gl0BchbAY/d264ALz0krp/JL4U+iI1bNkSun8e\nfRR27fp6LSBNApM4UOiLHIQ7vPFGCP9p06Bbt9D6v/xyLQctpUuhL1ILe/bAvHnhApBOw8UXhwtA\n//7q/pHSEnnom1kr4BGgO+DAOGAt8ARwIrAeGOXuO2ocp9CXSGzdGlr+FRWhK6h6+YeBA6F9+6hL\nJ3JoxRD6FcDz7j7ZzBoCRwP/Cnzs7v9hZj8CWrv7hBrHKfQlcuvWfb0ExOLFcOyxIfxHjgy/ArQa\nqBSbSEPfzFoCy939pBqvrwb6ufsWM2sLpN391BqfUehLUamqgjffDMtBz5gBf/4zjB4d1gA66yw4\n8sioSygSfej3Ah4CVgE9gdeB24CN7t468xkDtlc/zzpWoS9F7a234Ne/hoUL4Y9/DMM/q7uBTjtN\n9wIkGlGH/t8DS4Bz3f1VM7sP2AX8U3bIm9l2dz+mxrF+1113ffU8lUqRSqXqXBaRQvr449D9U90V\n9OmnYU+AgQNh2DBo1y7qEkpcpdNp0un0V8/vvvvuSEO/LbDE3TtnnvcF7gROAvq7+2YzKwMWq3tH\n4mTDhnABeO65sDfwWWeFuQAjR0LTplGXTuKsGG7kvgBc7+5rzGwiUP1Pfpu7/7uZTQBa6UauxNUX\nX8CsWWE00CuvwKWXhglhffuqC0jyrxhCvydhyGYj4D3CkM0GwJNARzRkUxLko4/CbOCKirAY3OjR\n4e/kk6MumcRF5KFf5xMr9CXG3GH58rAfwLRpcMopofU/ahS0ahV16aSUKfRFitzevTB/frgALFwY\nNoEZMwaGDAmbw4gcDoW+SAnZvh2eeCJ0/6xfH7aDHDMGevWKumRSKhT6IiVqzZrQ+p86NXT5jBkT\nLgJlZVGXTIqZQl+kxFVVwQsvhAvAzJnQp0/o/7/kEmjSJOrSSbFR6IvEyGefheGfU6fCH/4A554b\nJoBdcEGYBSyi0BeJqe3b4fnnwwSw2bPDYnDVXUBt2kRdOomKQl8kAaqqwj4AFRUwZ074BTBmTPgF\n0KJF1KWT+qTQF0mYzz4Lff+//jW8/DJ07x66gK66Sl1ASaDQF0mwv/wFliwJ6/9Mnw7HHRdmAKsL\nKL4U+iICwJdfhi6gxx7bvwto+HCNAooThb6I/I3PPoNnngmjgF57DS67LFwA+vYFq3NcSDFQ6IvI\nIW3aFBaBe+wxLQIXBwp9EakVd3jjjTACaPr0EPrXXBPWAurcWb8ASoVCX0QO2969YfG3GTPCPIBG\njcJOYNW7gWkpiOKl0BeRnLjDu+9+vRNYOh1a/mPGhGGgxx8fdQklm0JfRPLqyy/DfsDVo4DOOy9c\nAC66CBo3jrp0otAXkYLZvRuefjrcB3jzzbAJTHk5nHOO7gFERaEvIvViw4YwC7iiInQJDRkS+v/7\n99duYPVJoS8i9ap6FNCiReEewJIl4QJQXq7dwOqDQl9EIvXJJ2E3sEcfDb8GRo+GsWOhW7eoSxZP\nCn0RKRrvvBO6fx57DNq3h3Hjwgggdf/kj0JfRIrOvn1QWQlTpoT5ABdeCOPHh/7/I46IunSlTaEv\nIkVt2zaYNg0mTYKdO0Prf9w46NAh6pKVJoW+iJSMZctC+M+YAT17hk1ghgyBHj00BLS2FPoiUnI+\n/zyM/lmwAObPD/sCDB0a/i68EJo2jbqExUuhLyIlb926sBHM3LlhKehrr4V/+AeNADoQhb6IxMqG\nDfDIIzB5MjRvHlr+I0ZoL4BqCn0RiaWqKli+HH73u7AUdMOGcMst4VdAkrt/iiL0zawB8Bqw0d0v\nNrNjgCeAE4H1wCh331HjGIW+iNSKe5j9+1//FWYA33AD3HwztGsXdcnqX66hn68Rs7cCq4DqFJ8A\nVLp7F+C5zHMRkToxg0GDwqqfS5bArl1hxM+YMfD221GXrrTkHPpm1h4YBjwCVF99hgMVmccVwIhc\nzyMiAmHHr/vvh/feg65dYfDg0O//zDPwxRdRl6745aOl/wvgDqAq67U27r4l83gL0CYP5xER+Urr\n1nDnnfD++2Hj91/9Ck44Iaz98+yzYVaw/K2c1sMzs4uAre6+3MxSB/qMu7uZHbDzfuLEiV89TqVS\npFIH/AoRkYNq3Dgs8TB+PGzeDL/5Ddx9d1j0bfRo+Jd/Ke3dv9LpNOl0Om/fl9ONXDO7BxgN7AMa\nAy2AZ4CzgJS7bzazMmCxu59a41jdyBWRglm3Ltz4ffxxuOkmuP12aNky6lLlLtIbue7+Y3fv4O6d\ngSuB37v7aGAOUJ75WDkwK5fziIgcrpNPDqH/+uvw4Yfw7W/D//yPun3yvd5dddP9p8BgM1sDDMg8\nFxGpd506hbX+n302rPvfqxf8/vdRlyo6mpwlIonhDrNnw223hRm+P/85tG0bdakOT7GM0xcRKXpm\nYUmHlSvDJi89eoTlHpLU/lRLX0QSa8WKMMqnTRt4+OFwISh2aumLiNTR6afDK6/AueeG9f3/7d9g\nx45vPq6UKfRFJNGOPDKE/WuvwaZNcMop8J//Gd8uH3XviIhkWbMGysuhrCyM+mnRIuoS7U/dOyIi\nedSlC6TTYRbvOefAu+9GXaL8UuiLiNRw1FFhItcPfwj9+sHLL0ddovxR946IyCEsWBDW8HnoIRg5\nMurS5N69k9OCayIicTdkSNi8/eKLYft2uO66qEuUG4W+iMg36N0bFi+GgQPDBK/x46MuUd0p9EVE\naqFLl7Bl48CB4XmpBr9CX0SklqqDf8AAaN4cvve9qEt0+HQjV0TkML3xBpx/ftiisW/f+j23xumL\niNSzXr1g6lS4/PLSG8ev0BcRqYMhQ+Cee8Km7KW0Xo+6d0REcnDLLWFnrpkzw8ieQlP3johIhH7+\nc/jTn+Dee6MuSe2opS8ikqMNG+Dss+Hppwt/Y1ctfRGRiJ14IkyaBNdeC7t2RV2aQ1NLX0QkT8aP\nh8aN4cEHC3eOXFv6Cn0RkTzZsQNOOy0M5+zfvzDnUPeOiEiRaNUqrMZ5/fXw2WdRl+bA1NIXEcmz\nMWPCZus/+1n+v1vdOyIiRWbz5tDNs2RJ2HM3n9S9IyJSZNq2hTvuCH/FRqEvIlIAt94KK1aEVTmL\niUJfRKQAGjcOs3Vvuw327Yu6NF9T6IuIFMjIkXD88YW5oVtXOYW+mXUws8VmttLM3jazf868foyZ\nVZrZGjNbaGat8lNcEZHSYQZTpsAvfgGvvhp1aYJcW/p7gR+4e3egD3CzmXUFJgCV7t4FeC7zXEQk\ncTp2hF/9Cq6+Gnbvjro0eR6yaWazgAcyf/3cfYuZtQXS7n5qjc9qyKaIJMb48eAeWv65KJpx+mbW\nCXgeOA34wN1bZ143YHv186zPK/RFJDF274auXWH2bOjdu+7fk2vo52VjdDNrBjwN3OruuyxrJwF3\ndzM7YLpPnDjxq8epVIpUKpWP4oiIFJ1mzcIwznvvhccfr/1x6XSadDqdt3Lk3NI3syOBucCz7n5f\n5rXVQMrdN5tZGbBY3TsiknQ7d8JJJ8Hy5aGvvy4inZGb6bqZBKyqDvyMOUB55nE5MCuX84iIxEHL\nljB2LPzyl9GVIaeWvpn1BV4AVgDVX3QnsBR4EugIrAdGufuOGseqpS8iifPBB3DGGfDHP4aLwOEq\nmhu5h31ihb6IJNTVV4ebubfffvjHKvRFRErM66/DpZeG1n6DBod3rFbZFBEpMWeeGZZnWLiw/s+t\n0BcRicCNN8L//m/9n1fdOyIiEdi1KwzbXLUKyspqf5y6d0RESlDz5jBqVO7LMhwutfRFRCLy2mvw\nve/Be+/BEbVsgqulLyJSos48E1q3hkWL6u+cCn0RkYiYwQ03wKRJ9XhOde+IiERn27awHs/GjaGf\n/5uoe0dEpIR961tw3nlhyeX6oNAXEYnYVVfBtGn1cy5174iIRGz3bmjfHtauheOOO/Rn1b0jIlLi\nmjWDYcPgqacKfy6FvohIEbj66vrp4lH3johIEdizB044AZYtO/SuWureERGJgUaNYORIePrpwp5H\noS8iUiQuuADmzy/sOdS9IyJSJHbuDKN4tm6FJk0O/Bl174iIxETLltCzJ7zwQuHOodAXESkiQ4bA\nggWF+36FvohIEVHoi4gkyJlnwpYt8OGHhfl+hb6ISBFp0AAGDSrcpukKfRGRIlPILh4N2RQRKTKb\nNkGPHmHoZsOG+7+nIZsiIjHTrh2UlcGbb+b/uxX6IiJFqGdPePvt/H+vQl9EpAiddlqJhb6ZDTWz\n1Wa21sx+VKjziIjEUffusHJl/r+3IKFvZg2AB4ChQDfgKjPrWohziYjEUam19M8G1rn7enffC8wA\nLinQuUREYqdzZ9i2DT79NL/fW6jQbwdkzyfbmHlNRERq4YgjoGtXWLUqz9+b36/7igbgi4jkqHv3\n/HfxNPzmj9TJJqBD1vMOhNb+fiZOnPjV41QqRSqVKlBxRERKT/fuMH9+mo0b03n7zoLMyDWzhsC7\nwEDgI2ApcJW7v5P1Gc3IFRE5hHnz4L779l+HJ9cZuQVp6bv7PjP7J2AB0ACYlB34IiLyzQrRvaO1\nd0REipR72E1rwwZo3Tq8prV3RERiygy6dcvvJC2FvohIEcv3JC2FvohIEcv3cgwKfRGRIqbQFxFJ\nEHXviIgkSFkZ3HMPVFXl5/s0ZFNEpIRoyKaIiNSaQl9EJEEU+iIiCaLQFxFJEIW+iEiCKPRFRBJE\noS8ikiAKfRGRBFHoi4gkiEJfRCRBFPoiIgmi0BcRSRCFvohIgij0RUQSRKEvIpIgCn0RkQRR6IuI\nJIhCX0QkQRT6IiIJotAXEUkQhb6ISILUOfTN7Gdm9o6ZvWlmz5hZy6z37jSztWa22szOz09RRUQk\nV7m09BcC3d29J7AGuBPAzLoBVwDdgKHAg2aWuF8U6XQ66iIUlOpX2uJcvzjXLR/qHMbuXunuVZmn\nrwDtM48vAaa7+153Xw+sA87OqZQlKO7/8FS/0hbn+sW5bvmQrxb4eGBe5vEJwMas9zYC7fJ0HhER\nyUHDQ71pZpVA2wO89WN3/23mM/8K7HH3aYf4Kq97EUVEJF/Mve55bGZjgRuAge7+l8xrEwDc/aeZ\n5/OBu9z9lRrH6kIgIlIH7m51PbbOoW9mQ4F7gX7u/nHW692AaYR+/HbAIuBkz+XqIiIieXHI7p1v\ncD/QCKg0M4Al7n6Tu68ysyeBVcA+4CYFvohIccipe0dEREpLJOPnzWxoZuLWWjP7URRlyCcz62Bm\ni81spZm9bWb/nHn9GDOrNLM1ZrbQzFpFXda6MrMGZrbczKpv4Mepbq3M7KnMZMNVZnZOzOr3g8y/\ny7fMbJqZHVXK9TOzyWa2xczeynrtoPUptcmiB6lf3ibD1nvom1kD4AHCxK1uwFVm1rW+y5Fne4Ef\nuHt3oA9wc6ZOE4BKd+8CPJd5XqpuJXTZVf80jFPdfgnMc/euwOnAamJSPzNrB9wCnOnuPYAGwJWU\ndv2mEPIj2wHrU6KTRQ9Uv7xNho2i8mcD69x9vbvvBWYQJnSVLHff7O5vZB7vBt4h3MQeDlRkPlYB\njIimhLkxs/bAMOARoHrUQFzq1hL4rrtPBnD3fe6+k5jUL6Mh0NTMGgJNgY8o4fq5+4vAJzVePlh9\nSm6y6IHql8/JsFGEfjvgw6znsZq8ZWadgDMI/2PauPuWzFtbgDYRFStXvwDuAKqyXotL3ToDfzaz\nKWa2zMweNrOjiUn93H0TYZTdB4Sw3+HulcSkflkOVp84ThbNaTJsFKEf2zvHZtYMeBq41d13Zb+X\nGcFUcnU3s4uAre6+nK9b+fsp1bplNAR6Aw+6e2/gM2p0dZRy/cysNaEV3IkQEM3M7Nrsz5Ry/Q6k\nFvUp2brmYzJsFKG/CeiQ9bwD+1+pSpKZHUkI/KnuPivz8hYza5t5vwzYGlX5cnAuMNzM3gemAwPM\nbCrxqBuEf3sb3f3VzPOnCBeBzTGp3yDgfXff5u77gGeA7xCf+lU72L/HmnnTPvNayclMhh0GXJP1\n8mHXL4rQfw04xcw6mVkjwk2IORGUI28sTFSYBKxy9/uy3poDlGcelwOzah5b7Nz9x+7ewd07E24A\n/t7dRxODukG4HwN8aGZdMi8NAlYCvyUG9QM2AH3MrEnm3+kgwg35uNSv2sH+Pc4BrjSzRmbWGTgF\nWBpB+XKSmQx7B3BJ9eoHGYdfP3ev9z/gAuBdwk2HO6MoQ57r05fQ3/0GsDzzNxQ4hjAjeQ3h7nur\nqMuaYz37AXMyj2NTN6An8CrwJqEl3DJm9ZtIGFzwFuEm55GlXD/CL86PgD2E+4PjDlUf4MeZrFkN\nDIm6/HWo33hgLeECXp0vD9a1fpqcJSKSIMU+XlVERPJIoS8ikiAKfRGRBFHoi4gkiEJfRCRBFPoi\nIgmi0BcRSRCFvohIgvx/65k7D0LFXhgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10de57c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(z - sea_level)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}