mcflugen/rafem-bmi-example.ipynb

## rafem-bmi-example.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Before we begin...\n",
    "\n",
    "## Install the CSDMS tools\n",
    "\n",
    "Create a new folder to hold the Python installation, which will contain the CSDMS tools, and grab the Miniconda installer.\n",
    "```bash\n",
    "> mkdir py-csdms\n",
    "> cd py-csdms\n",
    "> curl http://repo.continuum.io/miniconda/Miniconda-latest-MacOSX-x86_64.sh -o miniconda.sh\n",
    "```\n",
    "There are also Miniconda installers for [Linux](http://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh). Unfortunately, we are yet to build our software stack on Windows.\n",
    "\n",
    "Install a fresh Python with the installer you just downloaded and add it to your `PATH`.\n",
    "```bash\n",
    "> bash miniconda.sh -f -b -p $(pwd)/conda\n",
    "> export PATH=$(pwd)/conda/bin:$PATH\n",
    "```\n",
    "\n",
    "Install the necessary CSDMS packages and then ipython notebook (and matplotlib so we can look at some pictures).\n",
    "```bash\n",
    "> conda install csdms-rafem csdms-cem coupling -c csdms\n",
    "> conda install ipython-notebook matplotlib\n",
    "```\n",
    "\n",
    "## Links\n",
    "* [Rafem source code](https://github.com/mcflugen/avulsion-bmi/tree/add-bmi-metadata): The BMI implementation is in `riverbmi.py`.\n",
    "* [Rafem description on CSDMS](http://csdms.colorado.edu/wiki/Model_help:RAFEM): Detailed information on the Rafem model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interacting with the Coastline Evolution Model BMI using Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting some stuff up for the example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some magic that allows us to view images within the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here I define a convenience function for plotting the water depth and making it look pretty. You don't need to worry too much about it's internals for this tutorial. It just saves us some typing later on."
   ]
  },
  {
   "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": [
    "### Create a model and query it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the `Rafem` class, and instantiate it. In Python, a model with a BMI will have no arguments for its constructor. Note that although the class has been instantiated, it's not yet ready to be run. We'll get to that later!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from cmt.components import Rafem\n",
    "raf = Rafem()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several values that can be set when setting up a *Rafem* simulation. These are initial values for parameters that can only be set before the model begins stepping through time. To see what these parameters (and their units) are and their default values, use the `defaults` attribute of the cem object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('delta_slope', (0.0001, '-')),\n",
       " ('upstream_elevation', (100.0, 'm')),\n",
       " ('number_of_columns', (500, '-')),\n",
       " ('rate_of_sea_level_rise', (0.025, 'm/yr')),\n",
       " ('initial_sea_level', (0.0, 'm')),\n",
       " ('row_spacing', (10, 'm')),\n",
       " ('number_of_rows', (1000, '-')),\n",
       " ('run_duration', (100.0, 'y')),\n",
       " ('column_spacing', (10, 'm')),\n",
       " ('rate_of_inlet_rise', (0.005, 'm/yr'))]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raf.defaults"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use the `setup` method to stage a cem simulation and set any parameters that we want to change from their default values. The first argument is a folder in which to stage the new simulation (if it already exists, it will be overwritten). New parameters are passed as keywords. Now that the model is staged, we can initialize it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "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)\n",
    "#          delta_slope=.0003, upstream_elevation=10., rate_of_sea_level_rise=.015)\n",
    "raf.initialize('_run_rafem/input.yaml')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now look at the current set of paramters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('rate_of_inlet_rise', 0.005),\n",
       " ('upstream_elevation', 100.0),\n",
       " ('number_of_columns', 100),\n",
       " ('rate_of_sea_level_rise', 0.015),\n",
       " ('initial_sea_level', 0.0),\n",
       " ('delta_slope', 0.0001),\n",
       " ('run_duration', 100.0),\n",
       " ('column_spacing', 5.0),\n",
       " ('row_spacing', 5.0),\n",
       " ('number_of_rows', 300)]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raf.parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As will all BMI objects, we can query the names of values that can set as the model updates through time.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('land_surface__elevation',\n",
       " 'channel_exit__x_coordinate',\n",
       " 'channel_exit__y_coordinate')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raf.get_input_var_names()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or the names of values that are provided by the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('channel_centerline__x_coordinate',\n",
       " 'channel_centerline__y_coordinate',\n",
       " 'channel_centerline__elevation',\n",
       " 'channel_water_sediment~bedload__volume_flow_rate',\n",
       " 'channel_exit__x_coordinate',\n",
       " 'channel_exit__y_coordinate',\n",
       " 'land_surface__elevation',\n",
       " 'sea_water_surface__elevation',\n",
       " 'avulsion_record')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raf.get_output_var_names()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data type: float64\n",
      "Units: m\n",
      "Grid id: 0\n",
      "Number of elements in grid: 30000\n",
      "Type of grid: uniform_rectilinear_grid\n"
     ]
    }
   ],
   "source": [
    "print \"Data type: %s\" % raf.get_var_type('land_surface__elevation')\n",
    "print \"Units: %s\" % raf.get_var_units('land_surface__elevation')\n",
    "print \"Grid id: %d\" % raf.get_var_grid('land_surface__elevation')\n",
    "print \"Number of elements in grid: %d\" % raf.get_grid_size(0)\n",
    "print \"Type of grid: %s\" % raf.get_grid_type(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set values of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "grid = raf.get_var_grid('land_surface__elevation')\n",
    "spacing = raf.get_grid_spacing(grid)\n",
    "shape = raf.get_grid_shape(grid)\n",
    "\n",
    "z0 = raf.get_value('land_surface__elevation').reshape(shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just for fun, we'll create a function that uplifts the land surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def uplift(z, t, sea_level=0.):\n",
    "    (y, x) = np.meshgrid(np.arange(z.shape[0]) * spacing[0],\n",
    "                         np.arange(z.shape[1]) * spacing[1], indexing='ij')\n",
    "\n",
    "    dz = np.exp(- t / 365. * 5.) * (np.cos(np.sqrt(x ** 2 + y ** 2) / 1e6 * np.pi * 20.))\n",
    "    on_land = z > sea_level\n",
    "    z[on_land] += dz[on_land]\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((300, 100), array([300, 100], dtype=int32))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z0.shape, shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "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/HP98xMLlzCMCRMgkGDAirFKqBAJYRBIlJr\niX1aAS8oSu3T+ihYLwhY2/B6Wgr4iNa2tPWCjSBoVEC8JyiDQQtITCAhBJIXiZhAhoRhGMBcZub8\nnj/W2mf2nDkzs2dyLnvI7/16zWv2XnuvvX/nzJz9O3uvvdeSmeGcc85lUWh0AM455yYPTxrOOecy\n86ThnHMuM08azjnnMvOk4ZxzLjNPGs455zKredKQtFnSg5JWSbovlrVJWi7pUUnLJLWm1r9M0gZJ\n6yWdWev4nHPOZVePMw0DOszsODM7MZZdCiw3s6OBn8V5JB0DnAscA5wFXCfJz4accy4n6nVAVtn8\n2cCSOL0EeHucXgTcbGZ9ZrYZ2AiciHPOuVyo15nGHZLul/TBWNZuZl1xugtoj9OHAVtSdbcAL6lD\njM455zJorsM+TjGzJyXNApZLWp9eaGYmabS+TLyfE+ecy4maJw0zezL+3i7pVsLlpi5Js81sm6Q5\nwFNx9a3A4anqc2NZyRgJxjnn3AjMrLypYNxUyw4LJe0HNJnZc5L2B5YBVwALgafN7GpJlwKtZnZp\nbAi/iZBYXgLcARxpqSAlmZlRHBgAoFgspqYH6O/rC9MDAxSLxdJ08ru0vDgwdFlqO/19fUO2GdYr\nlqaHLE/tp79vT2nd62/6Jhec947h6+zZMySu8jiG7rvya8sad3GgWPF9Scf0s988wOmvPZaB/sFt\nF5oKFJrC/1ahEH83ieaWQqksvTyZbm5J1ysMq5fMp7eZTJe2PWT54PZu+t7/cME584etk8wP3386\nrqYh+6y070r1kvnmKU3D3ouR3pd/+ucf8Zm/+5Mh21NzWN5vRn/q89Zvxq44n16W/j3S8v64jV3J\n/1X5tlPb2VUsDtYF/uuizbz+Xf+7tO5A0RhI1R0opqZjebosXZ5eVqnekOVWeXmlGAA23/rfHL7o\nfaX1h2x/jBiHbb/SOhm2Mdo2k/V33fFNWs44d/hrrPC6Sb+PSXnZe2v/eN6w/VSDpKokjVqfabQD\nt0pK9vUNM1sm6X5gqaQLgc3AOQBmtk7SUmAd4f/7Q+bd8DrnXG7UNGmY2SbgdRXKuwlnG5XqXAlc\nWcu4nHPOTYw/A1Ejx73m2EaHkMkRc9rHXikHXnvM4WOvlAMLFhzV6BDGdNhrTmh0CJnMeOVrGx1C\nJoUj/qDRIdSVJ40amSxJ4+VzZjc6hExee8xLGx1CJqctOLrRIYzpsNe8vtEhZHLQq4ZdpMilppdP\njs96tXjScM45l5knDeecc5l50nDOOZeZJw3nnHOZedJwzjmXWabnNCS1A6cQOhT8PbAWuN/MijWM\nzTnnXM6MmjQknU4Y6+IQ4DeEPqKmEboyP1LSt4HPmVlvrQN1zjnXeGOdabwV+KCZPV6+QFIL8Dbg\nTOA7NYjNOedczoyaNMzsk6Ms6wNurXpEzjnncitrm8bBwHuBeak6ZmYX1Sgu55xzOZT17qkfAS8D\nHgTuB1bGnzFJapK0StL343ybpOWSHpW0TFJrat3LJG2QtF7SmeN7Kc4552otay+3U83sYxPcx8WE\nrs4PjPOXAsvN7BpJn4rzyVga5wLHEMfSkHS036HlnHP5kfVM40ZJfyVpTjxTaJPUNlYlSXMJjelf\nAZLBP84GlsTpJYQ7sQAWATebWZ+ZbQY2EgZjcs45lxNZzzR2A9cAnwaSb/4GvHyMep8HPgnMSJW1\nm1lXnO4iDNQE4RmQe1LrbSGccTjnnMuJrEnjE4RhV3dk3bCktwFPmdkqSR2V1jEzG2PMbx+1zznn\nciRr0tgA7Bzntt8InC3prYQHAmdIugHokjTbzLZJmkN4YBBgK5AeaWduLBtm8eLFWBwbecGCBZw6\nf/44Q3POuRe3zs5OOjs7q77drEnj98BqSXcSLlXBGLfcmtnlwOUAkk4DPmFm50u6BngfcHX8fVus\ncjtwk6RrCZeljgLuq7TtxYsXUxwYAKBYLJamnXPOBR0dHXR0dJTmr7jiiqpsN2vSuC3+JJeLxPgv\nHSXrXwUslXQhsBk4B8DM1klaSrjTqh/4kJn55SnnnMuRrEljrZndny6Q9KdZd2JmdwF3xeluYOEI\n610JXJl1u8455+or6y23X5L0mmRG0juBz9QmJOecc3mV9UzjL4DvSHoXcCqhS5E31ywq55xzuZQp\naZjZY/Hs4jbgt8BbzOz3NY3MOedc7ow1nsaasqI2wiWteyWZmf1hzSJzzjmXO2OdaWRu7HbOOffi\nN1bS2GFmz4+2gqQDzey5KsbknHMup8a6e+p7kj4naYGk/ZNCSa+QdKGkZcBZtQ3ROedcXox1prGQ\n0EvtXwNvjD3b9gOPAD8E3mtm22obonPOubwYa7hXIySHH9YnHOecc3mW9eE+55xzzpOGc8657Dxp\nOOecyyxz0pB0qqT3x+lZko4YY/1pku6VtFrSWkmLY3mbpOWSHpW0TFJrqs5lkjZIWi/pzAm+Juec\nczWSKWnEA/4lwGWxaApw42h1zGwXcLqZvQ54HXCWpJOAS4HlZnY08LM4j6RjgHOBYwi38V4nyc+E\nnHMuR7IelP8MWAS8AGBmW4EDx6qU6p9qCtBCGFPjbGBJLF8CvD1OLwJuNrM+M9sMbAROzBifc865\nOsiaNHabWTGZST/oNxpJBUmrgS5gmZndB7SbWVdcpQtoj9OHAVtS1bcQRvBzzjmXE1mTxrcl/RfQ\nKumvCJeVvjJWJTMrxstTc4GTJB1bttwYfQRAH7nPOedyZMyu0SUJ+BbwKuA54GjgM2a2POtOzOzZ\nOL74W4AuSbPNbJukOcBTcbWtwOGpanNj2TCLFy/GiuHEZ8GCBZw6f37WUJxzbp/Q2dlJZ2dn1beb\ndRCmH5nZscCyrBuWNBPoN7MeSdMJgzZdBdwOvA+4Ov6+LVa5HbhJ0rWEy1JHAfdV2vbixYspDgwA\nUCwWS9POOeeCjo4OOjo6SvNXXHFFVbY7ZtIwM5O0UtKJsU0iqznAEklNhMtg3zKzH0m6B1gq6UJg\nM3BO3M86SUuBdYT+rT4UL18555zLiaxnGicD75H0W+IdVIR8MuIgTGa2Bji+Qnk3oSPESnWuBK7M\nGJNzzrk6y5o03hJ/J9/8VYNYnHPO5Vymu6ficxOthGcs/hQ4KJY555zbh2R9IvxiwhPgswjPVdwo\n6aJaBuaccy5/sl6e+kvgJDN7AUDSVcA9wBdrFZhz9dDWvh8Azz+7u8GRODc5jKdvp+II085NKq0z\np4Xfs6bT3fV7ert3MaNt2pB1kmTinBsq65nG14B7Jd1CaAR/O3B9zaLKoPeZ7sFnNQYG6O/rA+CA\n1tbh63Y/TX9fH60zZ9U1RpcPM9qmUiiI3u7dzGibSs+OXcw8bH96tu8srdPbvYvWmdPp7d4FQHfX\n72lr32/IOs65jEnDzK6VdBcwn3AH1QVmtqqmkY2iZ8d2ZhzcBgx9uK9YHKC7axsz2g6ht/vp0rIZ\nbYcAsOOJraX1kjL34tbWvh+93bspFETrrGn0dofLUN1dv6dQGHoTYM+OnbTOnE6hSfR27yoljuef\n3dOI0J3LpaxnGgCrgW2xjkl6qZk9XpuwRtc6c9aIT4HPaDuEnu1PMaPtEIqxq5Fk3bb22aUzkp4d\nT8VlRWa0tU0ojiR59Wx/qrSvGW1t9HY/7UmpRtIH/hHXiWcMrTPD5afmliYAerbvotA0+t3ivd1D\n1+nZvjMknme8zcM5yJg0JH0E+AdCP1Hpo/VrahHU3mqddeiYXYskB/XiwAA9O7YDcMBBwy9tjaS3\n+2laZ86iZ8d2WmcdCkB/3x56tm8vLa90qcxN3Iy2qfRs30XrrGk0tzSVzhR6dsQDe2ybSM4YenZU\n59JS7zO7aZ05drJybl+Q9Uzjo8ArzezpWgbTKK0zZ1EsFunu2jZiu0fPjnD20t21DRhMOjMObiud\nZQC0zppFcaA4pM5okrOSnh3bh1xKS5al22J6u58esu8XuxltU0uXhma0TS0dtJPLTYm29v1Kl5KS\nNonyM4a91bNjF60zp8VLV5483L4ra9J4HOitZSB50DpzVqlNJC0kk0Pp7tpWusSVpZPE9KWySpKE\nkU5W3V3bOOCg1nDpq+0QigMDwxJVcjYz0ctqeTejbfDg3NY+vTRd3gaRSM4oqp0oyiXJonXWNJ7v\nmVg7x5Y9od7cKVOqFpdz9TRq0pD08Tj5GNAp6QdA8mkxM7u2lsE1Qlv77NKBG6C3u5u29tkUB4q0\nzhr/3Vetsw6lZ/tTHNB68JDywXaXgbj9kISSs47k7CeJqVgcvEMsnM0MlBLLi0n6DqZCU4GeHWG6\nuSU/I//2dofLVWoOMe3s3YMOaA597Owefjf6lj176I99bybJYuPukIBmt7TUJWbnqmWsM40DCXdL\nPQ78jjBs6xTCbbdj9kAr6XDg68Chcf0vmdkXJbURxuh4GbGnWzPriXUuAz5AaDu5yMwyd8deLckZ\nQrE48Uby8u11d22jddahpbu6Rmt3Se4MG0vrzFnsePKJ0vRk1zortEOMdEaRJ8lZR6FJTG+bSt9z\nfRQKYvf0Arvi2GL9ZmCDiaLfrJQ8jpw6lX4zNu7e7WcdblIZNWmY2eLystjV+QFm9myG7fcBf2tm\nqyUdAKyUtBx4P7DczK6R9CngUuBSSccA5wLHEMbUuEPS0emhZuul1FBerM5YHa2zDqV725ND7uqq\nhuQspdJltclkRts0erbvrOnlpVrpe66vND11jzEV4uvQmAlw7pQpbPbE4SaRrH1P3SRpRhwbfA2w\nTtIlY9Uzs21mtjpOPw88TEgGZwNL4mpLCA8LAiwCbjazvtgh4kbgxHG8nlxrmz2ndttun01v99P0\ndj9duhssb5InsSs9fZ1cktoXzZs6tdTW4VzeZb1Q/Adm1ks4uP8YmAecP54dSZoHHAfcC7SbWVdc\n1EXoBBHgMGBLqtoWQpJxGcxoO4QZbYfQOnMWvd3dpbutGqWtfXqMa2rpSewZB08Nz1DMCsuSrjz2\ndXOnTBlX4tgc20TGm2zuXrGQnX0+0qWbuKxJo1lSCyFpfN/M+sjQppGIl6a+C1xsZs+ll8XR+Ubb\nlo/eNwEz2tqY0XYIzzy1jWdTZx59e8Z/gE4f/FtnTh+2PEkApfmZ05nRNpXurp2l5xtKt8vGh+R6\ntu8sXZJywdwpU9i8Zw9b+vpGXW9LXx/zpk5l465do9ZZtuIMAO68eyErfvlmOlcsZP6pdzA9Puzo\n3ERkveX2vwgN1g8Cv4hnDVnaNIjJ5rvADWaWjAfeJWm2mW2TNIfw0CDAVuDwVPW5sWyIxYsXY7Fd\nYMGCBZw6f37Gl7HvOfjQcOfVc888TaFJtEzZj4H+nUjhGvrU6cbO5wcvGSVtCoWC6H1md+ngP/Ow\n/QZvOY1dbSR3NCVPTSf1Q2N2XLZj14h3PvV275oUjd71NG/KlFID+bwK7Ryb9+xhbksL/cCR06ax\nq1gs1fneijcNWffMU+/gx794E6fPv6PUAN9fjxfhcqGzs5POzs6qb1cTGYZbkoAmMxv1fzCutwR4\n2sz+NlV+TSy7WtKlQKuZJQ3hNxHaMV4C3AEcmR4rXJKZWaq/qaF9TyW3pRYHBoZ1I5Lu2DBp4K60\nnfRzGIPrFUvTQ5an9tPft6e0bvk+Suvs2TMkrvI4hu678mvLGnc6jvT7AjDQHxqdB/qbmDrdeL5n\n8IG1QlOBQpOGPAWdThKFgoYkl2S6uaWQKi8Mq5fMJ4kiPV3a9pDlhaHbrlAvmR++/3RcTUPqVtp3\npXrJfPOUptJ0pddQ8X1JbS+5NTd991Qyvys5mKeWJb837t7N7JaWUkN5KWGYlQ7+u5L/q/Jtp7az\nq1gckjRW3vMW0gaKxkCq7kAxNR3L02Xp8vSySvWGLLfKyyvFUGl7Q+bHiHHYtiqtk2Ebo22zUtxj\nvW7S72NSXvbe2j+eN2w/1SAJM9vrb2nj6XuqJB7Es3xpOQV4D/CgpKSDw8uAq4Clki4k3nIbt7tO\n0lJgXdz+h2wiWc2NSZrCQH844Ox8vvLdXD07/EygkeZNmRLOOKZODYnDn+lwOTChpJGVmd3NyO0m\nC0eocyVwZc2Ccm4SSW7FnRsvQTnXaPl5zNY551zuZX1O4xxJM+L0ZyTdKun42obmnHMub7KeaXzG\nzHolzQfOAL4K/EftwnLOOZdHWZNG8jTQ24Avm9kPCH1QOeec24dkTRpbJX2J0C/UDyVNG0dd55xz\nLxJZD/znAD8Fzoy90R4MfLJmUTnnnMulrLfczgZ+aGa7JJ0O/CGDHQ4655zbR2Q907gF6Jd0JKFL\nkbmEJ7edc87tQ7ImjWLsMuR/Af9qZp8EatfPt3POuVzKmjT2SHoX8F7gB7HM+zRwzrl9TNak8QHg\nj4B/MrNNkl4O3FC7sJxzzuVRpqRhZg8BnwDWSjoW+J2ZXV3TyJxzzuVO1m5EOoBHgX8HrgM2SDot\nQ73rJXVJWpMqa5O0XNKjkpZJak0tu0zSBknrJZ057lfjnHOuprJenrqW8IzGAjNbAJwJfD5Dva8B\nZ5WVXQosN7OjgZ/FeeJYGucCx8Q610nyBwidcy5HMg/3amaPJDNm9igZnvEwsxXAM2XFZzP4jMcS\nwhCyAIuAm82sz8w2AxsJgzE555zLiawP962U9BXgRkDAu4H7J7jPdjPritNdQHucPgy4J7XeFsLo\nfc4553Iia9L4a+DDwEVxfgWhbWOvmJlJGm1kGR91xjnncmTMpCGpGXjAzF4FfK4K++ySNNvMtkma\nAzwVy7cCh6fWmxvLhlm8eDEWx0ZesGABp86fX4WwnHPuxaOzs5POzs6qbzdLu0S/pEckvczMfluF\nfd4OvA+4Ov6+LVV+k6RrCZeljgLuq7SBxYsXUxwIvbUXi8XStHPOuaCjo4OOjo7S/BVXXFGV7Wa9\nPNUGPCTpPuCFWGZmdvZolSTdDJwGzJT0O+DvgauApZIuBDYTetDFzNZJWgqsA/qBD5n5oMjOOZcn\nWZPGZ8rmMx3MzeydIyxaOML6VwJXZozJOedcnY2aNCQdRbjbqbOsfD7wZA3jcs45l0NjPafxBaC3\nQnlvXOacc24fMlbSaDezB8sLY9kRtQnJOedcXo2VNFpHWTatmoE455zLv7GSxv2S/qq8UNIHgZW1\nCck551xejXX31EeBWyW9m8EkcQIwFfizWgbmnHMuf0ZNGvGp7TcCpwPHEm61/YGZ/bwewTnnnMuX\nLE+EG/Dz+OOcc24f5uNVOOecy8yThnPOucw8aTjnnMvMk4ZzzrnMcpc0JJ0lab2kDZI+1eh4nHPO\nDcpV0pDUBPwbcBZwDPBOSa9ubFQTs2rN2kaHkMljT25rdAiZPLDu8UaHkMldv3i00SGM6Yk1Ex2p\nub6eXb+60SFkMvDY5PisV0uukgZwIrDRzDabWR/wTWBRg2OakMmSNDY92TX2SjnwwLrfNTqETH7x\niw2NDmFMT6yZHJ059D7yQKNDyKS46aFGh1BXeUsaLwHSR4ctscw551wO5C1p+Eh9zjmXY8rTiKqS\nTgYWm9lZcf4yoGhmV6fWyU/Azjk3iZiZ9nYbeUsazcAjwBnAE8B9wDvN7OGGBuaccw7IPkZ4XZhZ\nv6QPAz8FmoCvesJwzrn8yNWZhnPOuXzLW0P4qPLy4J+kwyXdKekhSWslXRTL2yQtl/SopGWSWlN1\nLotxr5d0Zp3jbZK0StL38xqnpFZJ35H0sKR1kk7KaZx/G//mayTdJGlqHuKUdL2kLklrUmXjjkvS\nCfG1bZD0L3WI8bPxb/6ApFskHdTIGEeKM7Xs45KKktryGqekj8T3dK2kdHtwdeI0s0nxQ7hctRGY\nB7QAq4FXNyiW2cDr4vQBhHaYVwPXAJfE8k8BV8XpY2K8LTH+jUChjvF+DPgGcHucz12cwBLgA3G6\nGTgob3ESbv9+DJga578FvC8PcQKnAscBa1Jl44kruepwH3BinP4RcFaNY3xz8p4AVzU6xpHijOWH\nAz8BNgFteYyTMPbRcqAlzs+qdpyT6UwjNw/+mdk2M1sdp58HHiYcUM4mHPyIv98epxcBN5tZn5lt\nJvzBTqxHrJLmAm8FvgIkd07kKs747fJUM7seQtuWmT2btzijZmC/eNPGfoQbNhoep5mtAJ4pKx5P\nXCdJmgMcaGb3xfW+nqpTkxjNbLmZFePsvcDcRsY4UpzRtcAlZWV5i/NvgH+Ox0jMbHu145xMSSOX\nD/5JmkfI9vcC7WaWPGLdBbTH6cMI8SbqGfvngU8CxVRZ3uI8Atgu6WuSfiPpy5L2z1ucZrYV+Bzw\nOCFZ9JjZ8rzFmTLeuMrLt1LfeD9A+KZLhVgaGqOkRcAWM3uwbFGu4gSOAhZIukdSp6TXVzvOyZQ0\nctdiL+kA4LvAxWb2XHqZhXO90WKu+euR9DbgKTNbxeBZxtAgchAn4dv78cB1ZnY88AJw6ZAgchCn\npIMJ397nET5sB0h6z5AgchBnxZ2OHVdDSfo0sMfMbmp0LOUk7QdcDvxDurhB4YylGTjYzE4mfFlc\nWu0dTKaksZVwTTFxOEMzZF1JaiEkjBvM7LZY3CVpdlw+B3gqlpfHPjeW1dobgbMlbQJuBt4k6YYc\nxrmF8C3u13H+O4Qksi1ncS4ENpnZ02bWD9wC/FEO40yM5++8JZbPLSuvebySLiBcQn13qjhPMb6C\n8EXhgfhZmguslNSesziJ+74FIH6eipJmVjPOyZQ07geOkjRP0hTgXOD2RgQiScBXgXVm9oXUotsJ\nDaPE37elys+TNEXSEYRTyPuoMTO73MwON7MjgPOAn5vZ+TmMcxvwO0lHx6KFwEPA9/MUJ/Bb4GRJ\n0+P/wEJgXQ7jTIzr7xz/Dr0Kd64JOD9VpyYknUX4RrzIzHaVxZ6LGM1sjZm1m9kR8bO0BTg+XvrL\nTZzRbcCbAOLnaYqZ7ahqnNVsza/1D/DHhDuVNgKXNTCO+YQ2gtXAqvhzFtAG3AE8CiwDWlN1Lo9x\nrwfe0oCYT2Pw7qncxQm8Fvg18ADhm9JBOY1zMeHGhzWExuWWPMRJOJN8AthDaPt7/0TiAk6Ir20j\n8MUax/gBYAMhGSefo+saGWNZnLuT97Js+WPEu6fyFmf8f7wh7ncl0FHtOP3hPuecc5lNpstTzjnn\nGsyThnPOucw8aTjnnMvMk4ZzzrnMPGk455zLzJOGc865zDxpuElD0ttjt9SvTJXNq9SFdR1i6VDs\nar4O+/qCpPlxenO6W+5xbudPJF1R3ejcvsaThptM3gmsiL8nNUmZPnuSDgFOMrO7Y9GEH6wysx8C\nfypp+kS34ZwnDTcpxM4hTwH+ktAlSqV1psWech+MveV2xPILFAb4+bHCgETpgWkulPSIpHtj77r/\nWmG7pykMYrUqbveAuOgASd+OA97cmFr/jLjeg5K+Gru9Sc4SrpK0EniHpDMl/UrSSklLY8++5f4c\n+HGFmKbH13OhpJcpDKzztfhabpS0UNLd8fW+IVW1E3jbqG+2c6PwpOEmi0XAj81sA/C0pOMrrPN/\ngAEz+0PC2cgSSVPjstcC5wCvAc6V9BJJhwF/B5xESEivpPI3+Y8DHzKz4whdyOyM5ccBFxMGuHm5\npDdKmgZ8DTgnxtFMGOOAuO0dZnYC8DPg08AZcX4lYbCscm+My9IOJPQl9A0z+yqhx9VXAP8PeFX8\nOc/M5gOfIHQfkbifMHiPcxPiScNNFu8kjJRH/F3pEtUpwI0AZvYIoU+jowkH65+Z2XNmtpvQyeA8\nwoBId5lZj4Vea79N5S6vfwl8XtJHCN1OD8Ty+8zsCQt98awmjAvySkJPuBvjOkuABaltJa/hZEKy\n+ZWkVcB7gZdW2PccYHtqXsD3gOvN7MZU+SYzeyjG8hAhKQGsja81sZ3QrbtzE9Lc6ACcG0ts+D0d\nOFaSEYb+NULvqMNWH2Ezu1PTA4T//fKzipHGHLla0g+APwF+Kekt49xmuuyF1PRyM3vXCPEmdgLp\nNggD7iZ03nlzqjwdS5HQKWAynf6cT2PwTMm5cfMzDTcZ/AXwdTObZ6F76pcCmySVX2ZZQRyTIXYL\n/VJCj56VkoERetU9TVKrwvCtf06Fy1OSXhG/xV8T64x0GcsIvTDPk/SKWHY+cFeFde8FTknWk7S/\npKMqrPcwcGRZ2d8Dz0j69wrrj+VoQo+mzk2IJw03GZwH3FpW9t1Ynh6R7jqgIOlBwhjy77MwVnLF\nUevM7AngSsIYF3cDm4DeCvu/WNIaSQ8QvsEnDdOVtrmb0EX1t2Mc/cB/lq9vYezmC4Cb43Z/RUhG\n5X4IdKR3EetfDEyXdNUIr89GmO6I23RuQrxrdLdPk7S/mb0QzzRuAb5qZt9rdFxpklYAbzOzZ/dy\nO+2ExvOF1YnM7Ys8abh9mqTPEkbgmwb81Mw+2uCQhpF0IrDTzPbqspKk1xPG4X6wOpG5fZEnDeec\nc5l5m4ZzzrnMPGk455zLzJOGc865zCbdw33x4S7nnHPjZGYjPfya2aRLGgCcvhiUeu0SFMrmYbAs\n+Z2uU9Dw+Up1s9RJ1iuvWym2keqUb7+8zmj7rfTaRypP6o/2flWIoUmiqWx+yO/C0PmkbKz5SnUn\nVKcasVaoM2oMZfs54eSflj5QzRLTCgWaY73mWJbMTysUhpU1S0wrmx9Sp2y+WcL6ixSLRnEgfJcq\nDhjFotHfVxwyX1peNPr3DFAs2rDlpTpxfkidvoEh80PWT20j2W5/XzFuo1i5Ttn6ybJQVhxWVqnO\n4Pqp1zJQpKm5hUJTEwCFQoFCU1NpvrmlhUKhiUJTIS4Py5pbWsJ8U1OpTqksrlMoFErrNLek9zF0\nefOUKYP7LcR1mgo0t0wZso3h+yik1h/cx5A6hbJYK8WV2kYSB4DSx4C94JennHPOZeZJwznnXGae\nNJxzzmXmScM551xmnjScc85l5knDOedcZp40nHPOZeZJwznnXGaeNJxzzmXmScM551xmnjScc85l\n5kljop7uL1+uAAAJEklEQVTa2OgIYMv6RkcAQP9jaxsdAj3rVzc6BADW//KxRofAL1ZsaHQIPPjw\n7xodAgCPPbGt0SGwak3jPx/V5EljojxplAxsavyH4tmcJI1Hfrmp0SGw4m5PGonHnvSkUW2eNJxz\nzmXmScM551xmMptcYxr5IEzOOTcx1RiEadIlDeecc43jl6ecc85l5knDOedcZrlNGpLeIekhSQOS\nji9bdpmkDZLWSzozVX6CpDVx2b/UIKbXSbpH0ipJv5b0hrFiqgVJH5H0sKS1kq5uRAypfX5cUlFS\nW73jkPTZ+D48IOkWSQfVO4bU/s6K+9og6VO13l/c5+GS7oyfk7WSLorlbZKWS3pU0jJJrXWKpyl+\nNr7fiDgktUr6TvyfWCfppEa8F5L+Nv491ki6SdLUWsch6XpJXZLWpMpG3OdefT7MLJc/wKuAo4E7\ngeNT5ccAq4EWYB6wkcG2mfuAE+P0j4CzqhzTMuAtcfqPgTtHialQo/fldGA50BLnZ9U7hlQshwM/\nATYBbQ14L96cbBu4CriqEe8F0BT3MS/uczXw6lq+93G/s4HXxekDgEeAVwPXAJfE8k8l70sd4vkY\n8A3g9jhf1ziAJcAH4nQzcFADYngJ8BgwNc5/C3hfreMATgWOA9akyiruc28/H7k90zCz9Wb2aIVF\ni4CbzazPzDYTXvBJkuYAB5rZfXG9rwNvr3JYRcI/IkArsHWUmE6s8r4TfwP8s5n1AZjZ9gbEkLgW\nuKSsrG5xmNlyMyvG2XuBufWOIToR2Ghmm+Pf5Zsxhpoys21mtjpOPw88TDhonU04gBJ/V/tzMIyk\nucBbga8AyR06dYsjnmWeambXA5hZv5k9W88YUpqB/SQ1A/sBT9Q6DjNbATxTVjzSPvfq85HbpDGK\nw4AtqfkthA9KefnWWF5NHwU+K+lx4LPAZWPEVAtHAQviZbJOSa9vQAxIWgRsMbMHyxbVNY6UDxDO\nLhsRw0uA9CPQ9XrNJZLmEb5p3gu0m1lXXNQFtNchhM8DnyR8sUrUM44jgO2SvibpN5K+LGn/OseA\nmW0FPgc8TkgWPWa2vN5xRCPtc68+H83ViW1iJC0nnGKXu9zMvl/veGDUmD4NLAQ+ama3SnoHcD3h\nEkklE76XeYwYmoGDzezk2KayFHh5tWPIEMdlQPpa6Gj3f9fivSj9j0j6NLDHzG6qRQwZNPS+dUkH\nAN8FLjaz56TBP4WZWa2fbZL0NuApM1slqaPSOnWIoxk4Hviwmf1a0heAS+scA5IOJnzDnwc8C3xb\n0nvqHUe5DPvMHE9Dk4aZjXTAHc1WwrX0xFxCptzK4OWJpHwr4zRaTJK+bmYXxdnvEE7FR4pp3PvO\nGMPfALfE9X4dG6FnVjuG0eKQdCzhm90D8QA1F1gp6aRqxzHW/4ikCwiXRc5IFVf9vRhD+f4OZ+g3\nuZqR1EJIGDeY2W2xuEvSbDPbFi/bPlXjMN4InC3prcA0YIakG+ocxxbCme+v4/x3CF9sttX5vVgI\nbDKzpwEk3QL8UQPigJHf/736fEyWy1Ppb7G3A+dJmiLpCMLlmvvMbBvQG++YEHA+cFuFbe2NJySd\nFqffBCRtLhVjqvK+E7fFfSPpaGCKme2oZwxmttbM2s3sCDM7gvCBPT6eCtctDklnES6JLDKzXalF\n9fx7ANwPHCVpnqQpwLkxhpqK/+dfBdaZ2RdSi24nNL4Sf1f7czCEmV1uZofH/4XzgJ+b2fn1jCN+\n/n8XPxMQDt4PAd+vVwzRb4GTJU2Pf5+FwLoGxAEjv/979/moZgt+NX+APyNcJ94JbAN+nFp2OaHx\nZj3xbqZYfgKwJi77Yg1iOoVwgFgN/A9w3Fgx1SCGFuCG+DpXAh31jqFCTI8R756q83uxgfAhXRV/\nrmvUe0G4m+6RuM/L6vS+zye0IaxOvQdnAW3AHYQvNcuA1jr+L5zG4N1TdY0DeC3wa+ABwtn4QY14\nL4DFhJsS1hAaoFtqHQdwM6ENZU88br5/tH3uzefDuxFxzjmX2WS5POWccy4HPGk455zLzJOGc865\nzDxpOOecy8yThnPOucw8aTjnnMvMk4ZzzrnMPGm4Fw2FsVdWpX4uieWdkk6o8r4+Kml6av6HkmZU\ncdvnj2P9KZLuktRUjf07Nxp/uM+9aEh6zswOrFB+J/BxM/tNFfe1CXi9xT6GqrjdZsKT/sfZYLfv\nWer9PaF79tE6bXRur/mZhtunSDpT0q8krZS0VNL+CqPuLU2t06HB0ef+Q2GUxrWSFseyiwjdS98p\n6WexbLPi6IWSPqYwatsaSRfHsnkKI8p9KW7rp5KmVQjxTcBvkoQRz5KujTGsk/R6hVEKH5X0f1P1\nbgPeXfU3zLkynjTci8n0sstT70gvjL0Bfxo4w8xOIHyj/xhhJMSTUpebziX05QOhC/Y3EPo1Ok3S\nsWb2RUI/Px1mlvSua3EfJwAXEAa1ORn4oKTXxXWOBP7NzI4FeoA/r/Aakv7NEgbsjjH8J/A9wkBc\nxwIXxK64IXTO9wacq7GGdo3uXJXtNLPjRlgmwkH8GOBXsUv3KcCvzGxA0k8I3Xt/l9DV+idivXMl\nfZDwWZkT668dZR/zgVvMbCeUusY+ldCz6CYbHLRqJWHMhXKzCb2ipiU95q4FHrI4sI6kx4CXAs/E\n17BH0v5m9sII8Tm31zxpuH3NcjN7V4XybwIfBrqB+83shdht9McJbRfPSvoaYbyI0RhDu/IXgwPc\n7E6VDwDTGW5nhX0k9Ypl2ygSxidPTAXSXcQ7V3V+ecrtKwy4BzhF0isAYnvGUXH5XYSR3z7I4KWp\nGcALhHFa2gndnyeei8vL97ECeHscT2F/wrjMKxh9ZMO0hwmXscZF0iHADjMbGG9d58bDk4Z7MSlv\n07gyvdDCYFUXADdLegD4FfDKuKwI/IAwHsUPYtkDhDEq1gPfAO5Obe5LwE+ShvDUPlYB/00Y1OYe\n4MtxOzB8SM1Kty7+GFgwwuuzEeoAnJ7E7Vwt+S23zuVMbAe5xMw2jqPOd4FPjaeOcxPhZxrO5c+l\nhEb3TBTGCb/NE4arBz/TcM45l5mfaTjnnMvMk4ZzzrnMPGk455zLzJOGc865zDxpOOecy+z/AyFI\nEGfre5n4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107ec1d10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_coast(raf.get_grid_spacing(grid), z0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sea_level_grid = raf.get_var_grid('sea_water_surface__elevation')\n",
    "raf.get_grid_size(sea_level_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raf.get_value('sea_water_surface__elevation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Advance the model through time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Advance the model in time while while getting elevations every time step, passing them to the `uplift` function and setting them in *Rafem* before the next update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for _ in xrange(3000 * 365 / 5):\n",
    "    raf.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAADLCAYAAACS0c8UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYXFWZ7/Hvr6q7k0ATmhZOuAQmEYjIMDNcFFBujUZg\nHIbLmRG8oSjDPDMeBa9cjw488wwCHtFxZpgZFTgRhDFyk4tKgtIMiBCICSSEQHJI1AQTrkkI5NJd\n9Z4/1trVu6urq3d36tbk/TxPP7332re3qrv2W2vtvdeSmeGcc85lkWt2AM4558YPTxrOOecy86Th\nnHMuM08azjnnMvOk4ZxzLjNPGs455zKre9KQtFLSU5IWSJoXy7olzZX0nKQ5krpS618saZmkpZJO\nqHd8zjnnsmtETcOAHjM7xMwOj2UXAXPNbAbwiziPpAOBM4EDgZOAayV5bcg551pEo07IKps/BZgV\np2cBp8XpU4FbzKzPzFYCy4HDcc451xIaVdO4X9ITks6NZVPMbG2cXgtMidN7AqtS264C9mpAjM45\n5zJoa8AxjjKzP0jaDZgraWl6oZmZpGp9mXg/J8451yLqnjTM7A/x90uS7iA0N62VtLuZrZG0B/Bi\nXH01sHdq86mxrGSEBOOcc24YZlZ+qWDUVM8OCyXtAOTN7HVJOwJzgMuBmcArZnaVpIuALjO7KF4I\nv5mQWPYC7gf2s1SQkmzW1y8HoFgoUCwWS9MA/X19sbxAsRCXFQv09/UNbBPXDcuKpflCfx+5fI5c\nXuRy4b3N5UVbe27QfC4ncvkw39aeK03ncgPTd/5yPn994ruH7CuXC/sbvK/coGOFfSm134HpocdM\nx5kfum3Zesm2yfy3/+1+vnjeBwaWd+SH7Dc9XdpHbiCeXF5DjqO28Dr6zehP/Y/1m7E5zqeXpX9v\nrlD2vX/6OZ+65ET64342x797+Xr9ZWWbi8WB5XHZ/EdPBKBQNApxWaEYf5fNl5eVr1e+n2W3X8/b\nT//UoLJq+0+XD1knXTbMepXiGfLa0vGbwQOz4bgPDaycbJt6zaT3Vxxmunyd8nNJteXFYdZNl8/7\nCbzrlJFjHCnWkY473PKR3pdk2TNz4IAPjG0fRRsSk913KfUgqSZJo941jSnAHZKSY/3QzOZIegKY\nLekcYCVwBoCZLZE0G1hC+Ix/xrwbXuecaxl1TRpmtgI4uEL5q4TaRqVtrgCuqGdczjnnxsafgaiT\nA6bv0ewQMnnPEW9vdgiZHHrMfs0OIZNdDhjyHan1TPvjZkeQzZ7vaHYE2ey6b7MjaChPGnVywPQ9\nmx1CJu85Ynz8wx927DhJGu88pNkhjGy6J42a2m18fIZqxZOGc865zDxpOOecy8yThnPOucw8aTjn\nnMvMk4ZzzrnMMj2nIWkKcBShQ8E3gcXAE2ZWrLqhc865t5SqSUPS8YSxLt4G/IbQR9REQlfm+0n6\nMfBNM9tQ70Cdc84130g1jQ8C55rZ78oXSGoHTgZOAG6tQ2zOOedaTNWkYWZfqbKsD7ij5hE555xr\nWVmvaewCfAKYltrGzOy8OsXlnHOuBWW9e+qnwB8BTwFPAPPjz4gk5SUtkHR3nO+WNFfSc5LmSOpK\nrXuxpGWSlko6YXQvxTnnXL1l7eV2gpl9cYzHOJ/Q1flOcf4iYK6ZXS3pwjifjKVxJnAgcSwNSTP8\nDi3nnGsdWWsaN0n6W0l7xJpCt6TukTaSNJVwMf37QDL4xynArDg9i3AnFsCpwC1m1mdmK4HlhMGY\nnHPOtYisNY0twNXApUDyzd+AkfrV/hbwFWByqmyKma2N02sJAzVBeAbk0dR6qwg1Dueccy0ia9L4\nMmHY1Zez7ljSycCLZrZAUk+ldczMRhjz20ftc865FpI1aSwDNo1y3+8FTpH0QcIDgZMl3QislbS7\nma2RtAfhgUGA1cDeqe2nxrIh7rj/AQCsWGTGtH04YPq0UYbmnHNvbb29vfT29tZ8v1mTxpvAQkkP\nEJqqYIRbbs3sEuASAEnHAV82s7MkXQ18Ergq/r4zbnIXcLOkawjNUvsD8yrt+/SZxwNQLBQoFv06\nuXPOlevp6aGnp6c0f/nll9dkv1mTxp3xJ2kuEqNvOkrWvxKYLekcYCVwBoCZLZE0m3CnVT/wGTPz\n5innnGshWZPGYjN7Il0g6S+zHsTMHgQejNOvAjOHWe8K4Iqs+3XOOddYWW+5/a6kP0lmJH0E+Gp9\nQnLOOdeqstY0/hq4VdJHgWMIXYp8oG5ROeeca0mZkoaZPR9rF3cCvwVONLM36xqZc865ljPSeBqL\nyoq6CU1aj0kyM/vTukXmnHOu5YxU08h8sds559xb30hJ42Uz21htBUk7mdnrNYzJOedcixrp7qmf\nSPqmpGMl7ZgUStpX0jmS5gAn1TdE55xzrWKkmsZMQi+1fwe8N/Zs2w88C9wLfMLM1tQ3ROecc61i\npOFejZAc7m1MOM4551pZ1of7nHPOOU8azjnnsvOk4ZxzLrPMSUPSMZI+Fad3kzR9hPUnSnpM0kJJ\niyVdFsu7Jc2V9JykOZK6UttcLGmZpKWSThjja3LOOVcnmZJGPOFfAFwcizqAm6ptY2abgePN7GDg\nYOAkSUcAFwFzzWwG8Is4j6QDgTOBAwm38V4ryWtCzjnXQrKelE8HTgXeADCz1cBOI22U6p+qA2gn\njKlxCjArls8CTovTpwK3mFmfma0ElgOHZ4zPOedcA2RNGlvMrDREXvpBv2ok5SQtBNYCc8xsHjDF\nzNbGVdYCU+L0nsCq1OarCCP4OeecaxFZk8aPJf0n0CXpbwnNSt8faSMzK8bmqanAEZIOKltuVB8B\n0Efuc865FjJi1+iSBPwIOAB4HZgBfNXM5mY9iJmtj+OLnwislbS7ma2RtAfwYlxtNbB3arOpsWyI\nO+5/IOy3WGTGtH04YPq0rKE4V1VeouCjDLu3gN7eXnp7e2u+36yDMP3UzA4C5mTdsaRdgX4zWydp\nEmHQpiuBu4BPAlfF33fGTe4CbpZ0DaFZan9gXqV9nz7zeACKhQLFYrHSKs6NSUc+x6b+QrPDcG6b\n9fT00NPTU5q//PLLa7LfEZOGmZmk+ZIOj9ckstoDmCUpT2gG+5GZ/VTSo8BsSecAK4Ez4nGWSJoN\nLCH0b/WZ2HzlXEN05HOl35uKIXHkc6JQ8H9D5xJZaxpHAh+X9FviHVSEfDLsIExmtgg4tEL5q4SO\nECttcwVwRcaYnNsm+Zw4/Mg5/PpXH2BSe55C0UpNU0kC8aYq5wbLmjROjL+TT5DqEItzDXXYkfcB\ngxNE+TWNQnFgOkksm/q8+cptvzLdPRWfm+giPGPxl8DOscy5tySvYThXWaaahqTzgXOB2wm1jJsk\nfc/MvlPP4Jyrh6SGATD/0VCJLhRDLSOR1Di8mcq5wbI+p/E3wBFm9jUz+yrhGse59QvLufqb/+iJ\npQSSz4WEkU4O6SQCg5uqnNtejaZvp+Iw086NC+XV6iRhJL+ThJFODgULF8cLRSslllqqxz6dq6es\nF8JvAB6TlDRPnQZcX7eoMlj30ksAFIsF+vv6wnShwOTutw1Zd/3LL5HL55jUObmhMbrW0l9l2dZC\nsdQklc+plDhKtY3U16t8TlAcuHCez2vEZzvS+6y0jCIUvAMENw5kShpmdo2kB4GjCXdQnW1mC+oa\nWRUbXn2Frt12G/RwX7EQPrQvv7Cayd1vY8Orr5bKdt51N/r7+tiy6XUAcjkhdTQneNdSkmsaeQ3U\nMob79p8+6edzIs/gRFD+YGBHPsfWglfK3VtL1poGwEJgTdzGJO1jZr+rT1jVVapNJLp2+x+se+lF\nJnd3U4wf2GJ8UGvCpJ1CoikUUC7UTvI5KI7xDsqu3Saxcd0WJndPBEIy2rh+K127TmTdy5vHtlNX\n1eZikYm56q2q/bGZaXOG3gKSpql5j5xQOsEPaZ5Kah0VagtJk9aktjwFMya15dnUX6Ajnxty8TxJ\nNGEmHCedVDzJuPEg691TnwP+gdBPVPoU+yf1CGpbTe5+WylRDMeK7QAUi0UmdYYP98Z1oznGRNa9\ntInuKTuw4dWQIHK5HN1TJgHQtetENq7fOobo3XA2xoSxMdYgk+TQJrHZjLZY1iaxsVBgYi6XKXEA\nw56s89KQK39J8ii/42rQNmS7cJ5OLukmMb9by7WqrDWNzwPvMLNX6hlMs2zaGE4Yk7uHP9FP7p7I\nhlc3M7l7Im3tuVKi2PDqZnL5gRNHuobRPWXSiImjs2sCb77ex+TuCQBsfnOg5b1z5w5yebFxXdhH\n164hIW1cv2W0L3FcStcqNpvRmcvRD3Tm82wuFmmLJ+fNxSITk8QRy0aqjUBomhp4CjyeuK1yIoDB\niaV0t1VxYLvy+fLaRpaHAsdS28jnRIe8zyzXGFmTxu+ADfUMpBVseHVzaHIqO9FP7p5QShgbXt1M\nLqdBiWL4/W2pWuOY3D2x1Ly14dUt5PIqrd/Z1VHabnL3BNra86x7eRMQkkcur7ds8khqBxNzOTYm\n01KpZlEuSRRtVL/YnVZ6PqPCPtOJY2uhOGSd8use6SfJk76qKtZAcpVv4U1qLcNJ94VVSXKsUg0n\n1RWKN3e5WquaNCR9KU4+D/RKugdIzoBmZtfUM7hmKD/Rd+0aT+g5DalVZNrfa1uYvMuEColoIAGl\nl21cv5XJ3RPYuG5r6Vjp6bBOSDCdO0/gzdffWk1g62KzEoSmps5cbthkMVZJwqhmUK2irLko3TyV\nLE9L1q90C28l6VpLOvHAwLWS4ZT3mVW+T+dqbaQ6/E5AJ6GmMZcwbGtnLB9xuFdJe0t6QNLTkhZL\nOi+Wd0uaK+k5SXMkdaW2uVjSMklLJZ0w1he2LTa8toXJ3ROY3D2hJhe0QyIIF8u7dptEZ9eEUvNW\nxeO/mq0GsXH9Fjq7JpT2Pd6tKxToyucbdrzRPqw3KFHkNKS2kd5nuraSrJtsX61mka41dFRoYksn\nqElt+Yrl5a9vuOY258aiak3DzC4rL4tdnXea2foM++8DvmBmCyV1AvMlzQU+Bcw1s6slXQhcBFwk\n6UDgTOBAwpga90uakR5qtlGSE/doaxbD7y80fa17aRNt7aN5prK6jetCnJO7Jw66HjLeJAmj1rWK\nWijvwLBQNDraKncvks+p9MxHuUFJpmgVp4csq1SDKKQ6WSy7DRggT4ghfQtwXqLAQLcoW8d6y6Db\n7mU6e0m6WdLkODb4ImCJpAtG2s7M1pjZwji9EXiGkAxOAWbF1WYRHhYEOBW4xcz6YoeIy4HDR/F6\nWlq12sW22rh+C507d4SfrtZ8BuXl/pDUkrufEo2uYdTCpr7CkOsFyXz5tQWoXqsZcn0kNZ+XSvtL\nEtWk9vyQdYckl2Sb/MBHPD3t3Fhl/S/6YzPbQDi5/wyYBpw1mgNJmgYcAjwGTDGztXHRWmBKnN4T\nWJXabBUhybgMNq7fGn7WbS0lkGZKksS6QoF1hQK7trWxrlCgM58vJY7NxeK4SxiVJAmj9KxHWcIo\nv8awtVAcknSq3bmVJIzh1s1LQ5q+0hfEsyYw50aSNWm0SWonJI27zawPsvd5EJumbgPON7PX08vi\n6HzV9uX/4WOQJJDy5NHZNWHU+9php/bS72Q6TRMHn/TfkLExlSS68vlSYkh+J4kjy62xra78JLy1\nWKzYbJSsm162qa/A1v6yGkqFvq+S8uTaSKV1k2sn+ZwG7uiKyztyA4nDL5K7bZH1ltv/JAzN+hTw\n37HWkOWaBjHZ3AbcaGbJeOBrJe1uZmsk7UF4aBBgNbB3avOpsWyQO+5/AAArFpkxbR8OmD4t48vY\n/pRu290lJIvkFt/krqsddurgzde3Dkkmbe053ny9jx0628PvncJvoJSEFL/99m3qp3/CwMl/x6JK\ny6pdp+hs0WsYtZLUJDryOTb1FQZfPC8MXTcvDal9pG/9LV2zSC6qp/rJSi9Pb1fej5Y/NLj96O3t\npbe3t+b71ViG4ZYkIG9mVa+8xvVmAa+Y2RdS5VfHsqskXQR0mVlyIfxmwnWMvYD7gf3SY4VLsllf\nDwOkV+p7qr+vL5YXBnUjku7UsJhqUy8Wi6X5Qn8fuXyOXF7k4gcwlxdt7blB8+nnNNrac6XpXC43\n6MJ5sqx82+RC+MC+coOOFfal1H5Vcb/JsoFt80O3LVsv2bZz5wnk8goJI3XrblvHwD46d+7gzY19\ng45R2kfqBJXLa8hxkqTRbzYoMfSbsTnOp5elf2+uUFZaP+4n/aR3peXprkRKy+Oy9Bgale56Ss+X\nl5WvN2Q/6XUrNAklJ/bkQb/0/qo1G6V72y3frlL3J+UxJuuUHgBMr5usV6msvLzaOuXnkmrLy19r\npRiyxjhSrCMdd7jlI70vWY6R5b0ti8nuu5R6kISZbXM1czR9T5XEk3iWW3WOAj4OPCUp6eDwYuBK\nYLakcwg1mDPifpdImg0sifv/jI0lq7kRJUkieUgwV6HJYuP6rTW7e8wFo33YbmtxoOkqaZ5KmrNg\nYPCoQRfFk1pIrLnkc/JahquZMSWNrMzsYYa/bjJzmG2uAK6oW1DOtbhKA0Glk0K6Saq0TXHoXVfp\ndZ2rlfF/FdK57UDp2kb/0P6v0ndupbs9KT2jkfNahqudrM9pnCFpcpz+qqQ7JB1a39Ccc+lbc7f2\nFzPXHMpvAa50kd25scha0/iqmW2QdDTwfuA64N/rF5ZzLq387igYPB6H1yRco2RNGsntRicD3zOz\newj9UDnnWpx3me5qKWvSWC3pu4R+oe6VNHEU2zrnxqDSLbvpmsVoahdbMw5G5dxIsp74zwDuA04w\ns3XALsBX6haVc64qvyPKNUvWW253B+41s82Sjgf+lIEOB51zDVBe23CuGbLWNG4H+iXtR+hSZCrh\nyW3nnHPbkaxJoxi7DPmfwL+Y2VeAPeoXlnPOuVaUNWlslfRR4BPAPbFsaHenzjnn3tKyJo1PA+8B\n/snMVkh6O3Bj/cJyzjnXijIlDTN7GvgysFjSQcDvzeyqukbmnHOu5WTtRqQHeA74N+BaYJmk4zJs\nd72ktZIWpcq6Jc2V9JykOZK6UssulrRM0lJJJ4z61TjnnKurrM1T1xCe0TjWzI4FTgC+lWG7G4CT\nysouAuaa2QzgF3GeOJbGmcCBcZtrJfkDhM4510IyD/dqZs8mM2b2HBme8TCzh4DXyopPYeAZj1mE\nIWQBTgVuMbM+M1sJLCcMxuScc65FZH24b76k7wM3AQI+BjwxxmNOMbO1cXotMCVO7wk8mlpvFWH0\nPueccy0ia9L4O+CzwHlx/iHCtY1tYmYmqVp/CN5XgnPOtZARk4akNuBJMzsA+GYNjrlW0u5mtkbS\nHsCLsXw1sHdqvamxbIg77n8AACsWmTFtHw6YPq0GYTnn3FtHb28vvb29Nd9vlusS/ZKelfRHZvbb\nGhzzLuCTwFXx952p8pslXUNoltofmFdpB6fPPB6AYqFA0XvvdM65IXp6eujp6SnNX3755TXZb9bm\nqW7gaUnzgDdimZnZKdU2knQLcBywq6TfA18DrgRmSzoHWEnoQRczWyJpNrAE6Ac+Y+YjyzjnXCvJ\nmjS+Wjaf6WRuZh8ZZtHMYda/ArgiY0zOOecarGrSkLQ/4W6n3rLyo4E/1DEu55xzLWik5zS+DWyo\nUL4hLnPOObcdGSlpTDGzp8oLY9n0+oTknHOuVY2UNLqqLJtYy0Ccc861vpGSxhOS/ra8UNK5wPz6\nhOScc65VjXT31OeBOyR9jIEkcRgwATi9noE555xrPVWTRnxq+73A8cBBhFtt7zGzXzYiOOecc60l\nyxPhBvwy/jjnnNuO+XgVzjnnMvOk4ZxzLjNPGs455zLzpOGccy6zlksakk6StFTSMkkXNjse55xz\nA1oqaUjKA/8KnAQcCHxE0jubG9XYLF3xQrNDyOTXj/2/ZoeQyfz/Xt7sEDJ57ZkFzQ5hZCuebnYE\n2bzwbLMjyOal8fEZqpWWShrA4cByM1tpZn3AfwGnNjmmMVm6Ynx0Avzrx55vdgiZ/OahcZI0li5s\ndggjW+lJo6Ze9qTRTHsBv0/Nr4plzjnnWkCrJQ0fqc8551qYWmlEVUlHApeZ2Ulx/mKgaGZXpdZp\nnYCdc24cMTNt6z5aLWm0Ac8C7wdeAOYBHzGzZ5oamHPOOSD7GOENYWb9kj4L3Afkges8YTjnXOto\nqZqGc8651tZqF8KrapUH/yTtLekBSU9LWizpvFjeLWmupOckzZHUldrm4hj3UkknNDjevKQFku5u\n1TgldUm6VdIzkpZIOqJF4/xC/JsvknSzpAmtEKek6yWtlbQoVTbquCQdFl/bMkn/3IAYvxH/5k9K\nul3Szs2Mcbg4U8u+JKkoqbtV45T0ufieLpaUvh5cmzjNbFz8EJqrlgPTgHZgIfDOJsWyO3BwnO4k\nXId5J3A1cEEsvxC4Mk4fGONtj/EvB3INjPeLwA+Bu+J8y8UJzAI+HafbgJ1bLU7C7d/PAxPi/I+A\nT7ZCnMAxwCHAolTZaOJKWh3mAYfH6Z8CJ9U5xg8k7wlwZbNjHC7OWL438HNgBdDdinESxj6aC7TH\n+d1qHed4qmm0zIN/ZrbGzBbG6Y3AM4QTyimEkx/x92lx+lTgFjPrM7OVhD/Y4Y2IVdJU4IPA94Hk\nzomWijN+uzzGzK6HcG3LzNa3WpxRG7BDvGljB8ING02P08weAl4rKx5NXEdI2gPYyczmxfV+kNqm\nLjGa2VwzK8bZx4CpzYxxuDija4ALyspaLc6/B74ez5GY2Uu1jnM8JY2WfPBP0jRCtn8MmGJma+Oi\ntcCUOL0nId5EI2P/FvAVoJgqa7U4pwMvSbpB0m8kfU/Sjq0Wp5mtBr4J/I6QLNaZ2dxWizNltHGV\nl6+msfF+mvBNlwqxNDVGSacCq8zsqbJFLRUnsD9wrKRHJfVKelet4xxPSaPlrthL6gRuA843s9fT\nyyzU9arFXPfXI+lk4EUzW8BALWNwEC0QJ+Hb+6HAtWZ2KPAGcNGgIFogTkm7EL69TyN82DolfXxQ\nEC0QZ8WDjhxXU0m6FNhqZjc3O5ZyknYALgH+IV3cpHBG0gbsYmZHEr4szq71AcZT0lhNaFNM7M3g\nDNlQktoJCeNGM7szFq+VtHtcvgfwYiwvj31qLKu39wKnSFoB3AK8T9KNLRjnKsK3uMfj/K2EJLKm\nxeKcCawws1fMrB+4HXhPC8aZGM3feVUsn1pWXvd4JZ1NaEL9WKq4lWLcl/BF4cn4WZoKzJc0pcXi\nJB77doD4eSpK2rWWcY6npPEEsL+kaZI6gDOBu5oRiCQB1wFLzOzbqUV3ES6MEn/fmSr/sKQOSdMJ\nVch51JmZXWJme5vZdODDwC/N7KwWjHMN8HtJM2LRTOBp4O5WihP4LXCkpEnxf2AmsKQF40yM6u8c\n/w4bFO5cE3BWapu6kHQS4RvxqWa2uSz2lojRzBaZ2RQzmx4/S6uAQ2PTX8vEGd0JvA8gfp46zOzl\nmsZZy6v59f4B/pxwp9Jy4OImxnE04RrBQmBB/DkJ6AbuB54D5gBdqW0uiXEvBU5sQszHMXD3VMvF\nCfwZ8DjwJOGb0s4tGudlhBsfFhEuLre3QpyEmuQLwFbCtb9PjSUu4LD42pYD36lzjJ8GlhGScfI5\nuraZMZbFuSV5L8uWP0+8e6rV4oz/jzfG484Hemodpz/c55xzLrPx1DzlnHOuyTxpOOecy8yThnPO\nucw8aTjnnMvMk4ZzzrnMPGk455zLzJOGGzcknRa7pX5HqmxapS6sGxBLj2JX8w041rclHR2nV6a7\n5R7lfv5C0uW1jc5tbzxpuPHkI8BD8fe4JinTZ0/S24AjzOzhWDTmB6vM7F7gLyVNGus+nPOk4caF\n2DnkUcDfELpEqbTOxNhT7lOxt9yeWH62wgA/P1MYkCg9MM05kp6V9FjsXfdfKuz3OIVBrBbE/XbG\nRZ2SfhwHvLkptf7743pPSboudnuT1BKulDQf+JCkEyQ9Imm+pNmxZ99yfwX8rEJMk+LrOUfSHykM\nrHNDfC03SZop6eH4et+d2rQXOLnqm+1cFZ403HhxKvAzM1sGvCLp0Arr/C+gYGZ/SqiNzJI0IS77\nM+AM4E+AMyXtJWlP4H8DRxAS0juo/E3+S8BnzOwQQhcym2L5IcD5hAFu3i7pvZImAjcAZ8Q42ghj\nHBD3/bKZHQb8ArgUeH+cn08YLKvce+OytJ0IfQn90MyuI/S4ui/wf4AD4s+Hzexo4MuE7iMSTxAG\n73FuTDxpuPHiI4SR8oi/KzVRHQXcBGBmzxL6NJpBOFn/wsxeN7MthE4GpxEGRHrQzNZZ6LX2x1Tu\n8vpXwLckfY7Q7XQhls8zsxcs9MWzkDAuyDsIPeEuj+vMAo5N7St5DUcSks0jkhYAnwD2qXDsPYCX\nUvMCfgJcb2Y3pcpXmNnTMZanCUkJYHF8rYmXCN26Ozcmbc0OwLmRxAu/xwMHSTLC0L9G6B11yOrD\n7GZLarpA+N8vr1UMN+bIVZLuAf4C+JWkE0e5z3TZG6npuWb20WHiTWwC0tcgDHiY0HnnLanydCxF\nQqeAyXT6cz6RgZqSc6PmNQ03Hvw18AMzm2ahe+p9gBWSyptZHiKOyRC7hd6H0KNnpWRghF51j5PU\npTB8619RoXlK0r7xW/zVcZvhmrGM0AvzNEn7xrKzgAcrrPsYcFSynqQdJe1fYb1ngP3Kyr4GvCbp\n3yqsP5IZhB5NnRsTTxpuPPgwcEdZ2W2xPD0i3bVATtJThDHkP2lhrOSKo9aZ2QvAFYQxLh4GVgAb\nKhz/fEmLJD1J+AafXJiutM8thC6qfxzj6Af+o3x9C2M3nw3cEvf7CCEZlbsX6EkfIm5/PjBJ0pXD\nvD4bZron7tO5MfGu0d12TdKOZvZGrGncDlxnZj9pdlxpkh4CTjaz9du4nymEi+czaxOZ2x550nDb\nNUnfIIzANxG4z8w+3+SQhpB0OLDJzLapWUnSuwjjcD9Vm8jc9siThnPOucz8moZzzrnMPGk455zL\nzJOGc865zMbdw33x4S7nnHOjZGbDPfya2bhLGgAcfxko9dolyJXNw0BZ8ju9TU5D5yttm2WbZL3y\nbSvFNtwYWLzlAAAKUklEQVQ25fsv36bacSu99uHKk+2rvV8VYshL5MvmB/3ODZ5Pykaar7TtmLap\nRawVtqkaQ9lxDjvyvtIHqk1iYi5HW9yuLZYl8xNzuSFlbRITy+YHbVM23yZh/UWKRaNYCN+ligWj\nWDT6+4qD5kvLi0b/1gLFog1ZXtomzg/apq8waH7Q+ql9JPvt7yvGfRQrb1O2frIslBWHlFXaZmD9\n1GspFMm3tZPL5wHI5XLk8vnSfFt7O7lcnlw+F5eHZW3t7WE+ny9tUyqL6+RyudI6be3pYwxe3tbR\nMXDcXFwnn6OtvWPQPoYeI5daf+AYg7bJlcVaKa7UPpI4AJQ+B2wDb55yzjmXmScN55xzmXnScM45\nl5knDeecc5l50nDOOZeZJw3nnHOZedJwzjmXmScN55xzmXnScM45l5knDeecc5l50nDOOZeZJ42x\nenF5syOAVUubHQEA/c8vbnYIrFu6sNkhALD0V883OwT++6FlzQ6Bp575fbNDAOD5F9Y0OwQWLGr+\n56OWPGmMlSeNksKK5n8o1rdI0nj2VyuaHQIPPexJI/H8Hzxp1JonDeecc5l50nDOOZeZzMbXmEY+\nCJNzzo1NLQZhGndJwznnXPN485RzzrnMPGk455zLrGWThqQPSXpaUkHSoWXLLpa0TNJSSSekyg+T\ntCgu++c6xHSwpEclLZD0uKR3jxRTPUj6nKRnJC2WdFUzYkgd80uSipK6Gx2HpG/E9+FJSbdL2rnR\nMaSOd1I81jJJF9b7ePGYe0t6IH5OFks6L5Z3S5or6TlJcyR1NSiefPxs3N2MOCR1Sbo1/k8skXRE\nM94LSV+If49Fkm6WNKHecUi6XtJaSYtSZcMec5s+H2bWkj/AAcAM4AHg0FT5gcBCoB2YBixn4NrM\nPODwOP1T4KQaxzQHODFO/znwQJWYcnV6X44H5gLtcX63RseQimVv4OfACqC7Ce/FB5J9A1cCVzbj\nvQDy8RjT4jEXAu+s53sfj7s7cHCc7gSeBd4JXA1cEMsvTN6XBsTzReCHwF1xvqFxALOAT8fpNmDn\nJsSwF/A8MCHO/wj4ZL3jAI4BDgEWpcoqHnNbPx8tW9Mws6Vm9lyFRacCt5hZn5mtJLzgIyTtAexk\nZvPiej8ATqtxWEXCPyJAF7C6SkyH1/jYib8Hvm5mfQBm9lITYkhcA1xQVtawOMxsrpkV4+xjwNRG\nxxAdDiw3s5Xx7/JfMYa6MrM1ZrYwTm8EniGctE4hnECJv2v9ORhC0lTgg8D3geQOnYbFEWuZx5jZ\n9QBm1m9m6xsZQ0obsIOkNmAH4IV6x2FmDwGvlRUPd8xt+ny0bNKoYk9gVWp+FeGDUl6+OpbX0ueB\nb0j6HfAN4OIRYqqH/YFjYzNZr6R3NSEGJJ0KrDKzp8oWNTSOlE8TapfNiGEvIP0IdKNec4mkaYRv\nmo8BU8xsbVy0FpjSgBC+BXyF8MUq0cg4pgMvSbpB0m8kfU/Sjg2OATNbDXwT+B0hWawzs7mNjiMa\n7pjb9Ploq01sYyNpLqGKXe4SM7u70fFA1ZguBWYCnzezOyR9CLie0ERSyZjvZR4hhjZgFzM7Ml5T\nmQ28vdYxZIjjYiDdFlrt/u96vBel/xFJlwJbzezmesSQQVPvW5fUCdwGnG9mr0sDfwozs3o/2yTp\nZOBFM1sgqafSOg2Iow04FPismT0u6dvARQ2OAUm7EL7hTwPWAz+W9PFGx1EuwzEzx9PUpGFmw51w\nq1lNaEtPTCVkytUMNE8k5asZpWoxSfqBmZ0XZ28lVMWHi2nUx84Yw98Dt8f1Ho8XoXetdQzV4pB0\nEOGb3ZPxBDUVmC/piFrHMdL/iKSzCc0i708V1/y9GEH58fZm8De5upHUTkgYN5rZnbF4raTdzWxN\nbLZ9sc5hvBc4RdIHgYnAZEk3NjiOVYSa7+Nx/lbCF5s1DX4vZgIrzOwVAEm3A+9pQhww/Pu/TZ+P\n8dI8lf4WexfwYUkdkqYTmmvmmdkaYEO8Y0LAWcCdFfa1LV6QdFycfh+QXHOpGFONj524Mx4bSTOA\nDjN7uZExmNliM5tiZtPNbDrhA3torAo3LA5JJxGaRE41s82pRY38ewA8AewvaZqkDuDMGENdxf/z\n64AlZvbt1KK7CBdfib9r/TkYxMwuMbO94//Ch4FfmtlZjYwjfv5/Hz8TEE7eTwN3NyqG6LfAkZIm\nxb/PTGBJE+KA4d//bft81PIKfi1/gNMJ7cSbgDXAz1LLLiFcvFlKvJsplh8GLIrLvlOHmI4inCAW\nAr8GDhkppjrE0A7cGF/nfKCn0TFUiOl54t1TDX4vlhE+pAviz7XNei8Id9M9G495cYPe96MJ1xAW\npt6Dk4Bu4H7Cl5o5QFcD/xeOY+DuqYbGAfwZ8DjwJKE2vnMz3gvgMsJNCYsIF6Db6x0HcAvhGsrW\neN78VLVjbsvnw7sRcc45l9l4aZ5yzjnXAjxpOOecy8yThnPOucw8aTjnnMvMk4ZzzrnMPGk455zL\nzJOGc865zDxpuLcMhbFXFqR+LojlvZIOq/GxPi9pUmr+XkmTa7jvs0axfoekByXla3F856rxh/vc\nW4ak181spwrlDwBfMrPf1PBYK4B3WexjqIb7bSM86X+IDXT7nmW7rxG6Z6/WaaNz28xrGm67IukE\nSY9Imi9ptqQdFUbdm51ap0cDo8/9u8IojYslXRbLziN0L/2ApF/EspWKoxdK+qLCqG2LJJ0fy6Yp\njCj33biv+yRNrBDi+4DfJAkj1pKuiTEskfQuhVEKn5P0j6nt7gQ+VvM3zLkynjTcW8mksuapD6UX\nxt6ALwXeb2aHEb7Rf5EwEuIRqeamMwl9+UDogv3dhH6NjpN0kJl9h9DPT4+ZJb3rWjzGYcDZhEFt\njgTOlXRwXGc/4F/N7CBgHfBXFV5D0r9ZwoAtMYb/AH5CGIjrIODs2BU3hM753o1zddbUrtGdq7FN\nZnbIMMtEOIkfCDwSu3TvAB4xs4KknxO6976N0NX6l+N2Z0o6l/BZ2SNuv7jKMY4GbjezTVDqGvsY\nQs+iK2xg0Kr5hDEXyu1O6BU1LekxdzHwtMWBdSQ9D+wDvBZfw1ZJO5rZG8PE59w286Thtjdzzeyj\nFcr/C/gs8CrwhJm9EbuN/hLh2sV6STcQxouoxhjclb8YGOBmS6q8AExiqE0VjpFsVyzbR5EwPnli\nApDuIt65mvPmKbe9MOBR4ChJ+wLE6xn7x+UPEkZ+O5eBpqnJwBuEcVqmELo/T7wel5cf4yHgtDie\nwo6EcZkfovrIhmnPEJqxRkXS24CXzaww2m2dGw1PGu6tpPyaxhXphRYGqzobuEXSk8AjwDvisiJw\nD2E8inti2ZOEMSqWAj8EHk7t7rvAz5ML4aljLAD+L2FQm0eB78X9wNAhNSvduvgz4NhhXp8Nsw3A\n8UncztWT33LrXIuJ10EuMLPlo9jmNuDC0Wzj3Fh4TcO51nMR4aJ7JgrjhN/pCcM1gtc0nHPOZeY1\nDeecc5l50nDOOZeZJw3nnHOZedJwzjmXmScN55xzmf1/C/ZAXwZxITYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10860cb10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "z = raf.get_value('land_surface__elevation').reshape(shape)\n",
    "plot_coast(raf.get_grid_spacing(0), z - raf.get_value('sea_water_surface__elevation') )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 45.])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raf.get_value('sea_water_surface__elevation')"
   ]
  }
 ],
 "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
}