diff --git a/lessons/pymt/cem.ipynb b/lessons/pymt/cem.ipynb
new file mode 100644
index 0000000..575fa35
--- /dev/null
+++ b/lessons/pymt/cem.ipynb
@@ -0,0 +1,480 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Coastline Evolution Model\n",
+ "\n",
+ "* Link to this notebook: https://github.com/csdms/pymt/blob/master/docs/demos/cem.ipynb\n",
+ "* Install command: `$ conda install notebook pymt_cem`\n",
+ "* Download local copy of notebook:\n",
+ "\n",
+ " `$ curl -O https://raw.githubusercontent.com/csdms/pymt/master/docs/demos/cem.ipynb`\n",
+ "\n",
+ "This example explores how to use a BMI implementation using the CEM model as an example.\n",
+ "\n",
+ "### Links\n",
+ "* [CEM source code](https://github.com/csdms/cem-old/tree/mcflugen/add-function-pointers): Look at the files that have *deltas* in their name.\n",
+ "* [CEM description on CSDMS](http://csdms.colorado.edu/wiki/Model_help:CEM): Detailed information on the CEM model."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Interacting with the Coastline Evolution Model BMI using Python"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Some magic that allows us to view images within the notebook."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%matplotlib inline\n",
+ "\n",
+ "from tqdm import tqdm"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Import the `Cem` 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": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pymt.models\n",
+ "cem = pymt.models.Cem()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Even though we can't run our waves model yet, we can still get some information about it. *Just don't try to run it.* Some things we can do with our model are get the names of the input variables."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.output_var_names"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.input_var_names"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can also get information about specific variables. Here we'll look at some info about wave direction. This is the main input of the Cem model. Notice that BMI components always use [CSDMS standard names](http://csdms.colorado.edu/wiki/CSDMS_Standard_Names). The CSDMS Standard Name for wave angle is,\n",
+ "\n",
+ " \"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\"\n",
+ "\n",
+ "Quite a mouthful, I know. With that name we can get information about that variable and the grid that it is on (it's actually not a one)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "angle_name = 'sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity'\n",
+ "\n",
+ "print(\"Data type: %s\" % cem.get_var_type(angle_name))\n",
+ "print(\"Units: %s\" % cem.get_var_units(angle_name))\n",
+ "print(\"Grid id: %d\" % cem.get_var_grid(angle_name))\n",
+ "print(\"Number of elements in grid: %d\" % cem.get_grid_number_of_nodes(0))\n",
+ "print(\"Type of grid: %s\" % cem.get_grid_type(0))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "OK. We're finally ready to run the model. Well not quite. First we initialize the model with the BMI **initialize** method. Normally we would pass it a string that represents the name of an input file. For this example we'll pass **None**, which tells Cem to use some defaults."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "args = cem.setup(number_of_rows=100, number_of_cols=200, grid_spacing=200.)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.initialize(*args)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Before running the model, let's set a couple input parameters. These two parameters represent the wave height and wave period of the incoming waves to the coastline."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "\n",
+ "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": "markdown",
+ "metadata": {},
+ "source": [
+ "The main output variable for this model is *water depth*. In this case, the CSDMS Standard Name is much shorter:\n",
+ "\n",
+ " \"sea_water__depth\"\n",
+ "\n",
+ "First we find out which of Cem's grids contains water depth. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid_id = cem.get_var_grid('sea_water__depth')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "With the *grid_id*, we can now get information about the grid. For instance, the number of dimension and the type of grid (structured, unstructured, etc.). This grid happens to be *uniform rectilinear*. If you were to look at the \"grid\" types for wave height and period, you would see that they aren't on grids at all but instead are scalars."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid_type = cem.get_grid_type(grid_id)\n",
+ "grid_rank = cem.get_grid_ndim(grid_id)\n",
+ "print('Type of grid: %s (%dD)' % (grid_type, grid_rank))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Because this grid is uniform rectilinear, it is described by a set of BMI methods that are only available for grids of this type. These methods include:\n",
+ "* get_grid_shape\n",
+ "* get_grid_spacing\n",
+ "* get_grid_origin"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "spacing = np.empty((grid_rank, ), dtype=float)\n",
+ "\n",
+ "shape = cem.get_grid_shape(grid_id)\n",
+ "cem.get_grid_spacing(grid_id, out=spacing)\n",
+ "\n",
+ "print('The grid has %d rows and %d columns' % (shape[0], shape[1]))\n",
+ "print('The spacing between rows is %f and between columns is %f' % (spacing[0], spacing[1]))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Allocate memory for the water depth grid and get the current values from `cem`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.empty(shape, dtype=float)\n",
+ "cem.get_value('sea_water__depth', out=z)"
+ ]
+ },
+ {
+ "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": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_coast(spacing, z):\n",
+ " import matplotlib.pyplot as plt\n",
+ " \n",
+ " xmin, xmax = 0., z.shape[1] * spacing[0] * 1e-3\n",
+ " ymin, ymax = 0., z.shape[0] * spacing[1] * 1e-3\n",
+ "\n",
+ " plt.imshow(z, extent=[xmin, xmax, ymin, ymax], origin='lower', cmap='ocean')\n",
+ " plt.colorbar().ax.set_ylabel('Water Depth (m)')\n",
+ " plt.xlabel('Along shore (km)')\n",
+ " plt.ylabel('Cross shore (km)')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "It generates plots that look like this. We begin with a flat delta (green) and a linear coastline (y = 3 km). The bathymetry drops off linearly to the top of the domain."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_coast(spacing, z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Right now we have waves coming in but no sediment entering the ocean. To add some discharge, we need to figure out where to put it. For now we'll put it on a cell that's next to the ocean."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Allocate memory for the sediment discharge array and set the discharge at the coastal cell to some value."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "qs = np.zeros_like(z)\n",
+ "qs[0, 100] = 1250"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The CSDMS Standard Name for this variable is:\n",
+ "\n",
+ " \"land_surface_water_sediment~bedload__mass_flow_rate\"\n",
+ "\n",
+ "You can get an idea of the units based on the quantity part of the name. \"mass_flow_rate\" indicates mass per time. You can double-check this with the BMI method function **get_var_units**."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.get_var_units('land_surface_water_sediment~bedload__mass_flow_rate')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.time_step, cem.time_units, cem.time"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Set the bedload flux and run the model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for time in tqdm(range(3000)):\n",
+ " cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
+ " cem.update_until(time)\n",
+ "\n",
+ "cem.get_value('sea_water__depth', out=z)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.time"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.get_value('sea_water__depth', out=z)\n",
+ "plot_coast(spacing, z)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "val = np.empty((5, ), dtype=float)\n",
+ "cem.get_value(\"basin_outlet~coastal_center__x_coordinate\", val)\n",
+ "val / 100."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's add another sediment source with a different flux and update the model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "qs[0, 150] = 1500\n",
+ "for time in tqdm(range(3750)):\n",
+ " cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
+ " cem.update_until(time)\n",
+ " \n",
+ "cem.get_value('sea_water__depth', out=z)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_coast(spacing, z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Here we shut off the sediment supply completely."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "qs.fill(0.)\n",
+ "for time in tqdm(range(4000)):\n",
+ " cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
+ " cem.update_until(time)\n",
+ " \n",
+ "cem.get_value('sea_water__depth', out=z)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_coast(spacing, z)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lessons/pymt/cem_and_waves.ipynb b/lessons/pymt/cem_and_waves.ipynb
new file mode 100644
index 0000000..a9f2b54
--- /dev/null
+++ b/lessons/pymt/cem_and_waves.ipynb
@@ -0,0 +1,602 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Couple two models in *pymt*\n",
+ "\n",
+ "In this exercise we will couple two component using *pymt*. We will couple the Coastline Evolution Model (CEM),\n",
+ "which models the transport of sediment alongshore, with a wave model that provides incoming wave direction and\n",
+ "characteristics.\n",
+ "\n",
+ "* Explore the base-case river simulation\n",
+ "* How does a river system respond to climate change?\n",
+ "* How do humans affect river sediment loads?\n",
+ "\n",
+ "## The Coastline Evolution Model (CEM)\n",
+ "\n",
+ "The Coastline Evolution Model (CEM) addresses predominately sandy, wave-dominated coastlines on time-scales ranging from years to millenia and on spatial scales ranging from kilometers to hundreds of kilometers. Shoreline evolution results from gradients in wave-driven alongshore sediment transport.\n",
+ "\n",
+ "At its most basic level, the model follows the standard *one-line* modeling approach, where the cross-shore dimension is collapsed into a single data point. However, the model allows the planview shoreline to take on arbitrary local orientations, and even fold back upon itself, as complex shapes such as capes and spits form under some wave climates (distributions of wave influences from different approach angles). The model works on a 2D grid.\n",
+ "\n",
+ "CEM has been used to represent varying geology underlying a sandy coastline and shoreface in a simplified manner and enables the simulation of coastline evolution when sediment supply from an eroding shoreface may be constrained. CEM also supports the simulation of human manipulations to coastline evolution through beach nourishment or hard structures.\n",
+ "\n",
+ "CEM authors & developers include:\n",
+ "* Andrew Ashton\n",
+ "* Brad Murray\n",
+ "* Jordan Slot\n",
+ "* Jaap Nienhuis and others.\n",
+ "\n",
+ "This version is adapted from a CSDMS teaching notebook, listed below. \n",
+ "It has been created by Irina Overeem, October 2019 for a Sedimentary Modeling course.\n",
+ "\n",
+ "### Key References\n",
+ "\n",
+ "Ashton, A.D., Murray, B., Arnault, O. 2001. Formation of coastline features by large-scale instabilities induced by high-angle waves, Nature 414.\n",
+ "\n",
+ "Ashton, A. D., and A. B. Murray (2006), High-angle wave instability and emergent shoreline shapes: 1. Modeling of sand waves, flying spits, and capes, J. Geophys. Res., 111, F04011, doi:10.1029/2005JF000422.\n",
+ "\n",
+ "\n",
+ "### Links\n",
+ "* [CEM source code](https://github.com/csdms/cem-old/tree/mcflugen/add-function-pointers): Look at the files that have *deltas* in their name.\n",
+ "* [CEM description on CSDMS](http://csdms.colorado.edu/wiki/Model_help:CEM): Detailed information on the CEM model."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Interacting with the Coastline Evolution Model BMI using Python"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from tqdm import tqdm"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Import the `Cem` model into your environment."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pymt.models\n",
+ "cem = pymt.models.Cem()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Even though we can't run our model yet, we can still get some information about it. Some things we can do with our model are to get help, and get the names of the input variables or output variables. In looking at the "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.input_var_names"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.output_var_names"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can also get information about specific variables. Here we'll look at some info about wave direction. This is the main input of the CEM model. What do you think the more conventional names for these variables are?\n",
+ "\n",
+ "| Conventional Name | Standard Name |\n",
+ "| :--------------------- | :------------------------------------------------------------------ |\n",
+ "| ??? | sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity |\n",
+ "| ??? | sea_surface_water_wave__period |\n",
+ "| ??? | sea_surface_water_wave__height |\n",
+ "\n",
+ "To help us out, we can get some additional information about each of the variables."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "angle_name = 'sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity'\n",
+ "\n",
+ "print(\"Data type: %s\" % cem.var_type(angle_name))\n",
+ "print(\"Units: %s\" % cem.var_units(angle_name))\n",
+ "print(\"Grid id: %d\" % cem.var_grid(angle_name))\n",
+ "print(\"Number of elements in grid: %d\" % cem.grid_node_count(0))\n",
+ "print(\"Type of grid: %s\" % cem.grid_type(0))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We now get the model ready for time stepping. Remember the lifecycle of the model is:\n",
+ "* *setup*\n",
+ "* *initialize*\n",
+ "* *update*\n",
+ "* *finalize*\n",
+ "\n",
+ "For this example we'll set up a simulation with a grid of 100 rows and 200 columns with a grid\n",
+ "resolution of 200.0."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "args = cem.setup(number_of_rows=100, number_of_cols=200, grid_spacing=200.)\n",
+ "cem.initialize(*args)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Changing a model's state\n",
+ "\n",
+ "Because the CEM model has input variables (unlike *HydroTrend* in the previous example), we\n",
+ "are able to change variables within the CEM model. The is done with the ***set_value***\n",
+ "method.\n",
+ "\n",
+ "For our first example we'll set the incoming wave height, period, and angle (in radians)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.set_value(\"sea_surface_water_wave__height\", 1.5)\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.0 * np.pi / 180.)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Grids\n",
+ "\n",
+ "This example is also different from the previous example in that it generates output that\n",
+ "is on a grid (as opposed to scalar data). The main output for CEM is *sea_water__depth*\n",
+ "and operates on a grid whose size we set when we *setup* this simulation.\n",
+ "\n",
+ "*pymt* models can have multiple grids. This allows for models to calculate some of\n",
+ "its state variables on scalars and others of 2D grids, for example. Models could also\n",
+ "maintain several grids of differing resolution. In *pymt* each grid has an ID\n",
+ "associated with it.\n",
+ "\n",
+ "We use the ***var_grid*** function to get the grid ID."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.var_grid(\"sea_water__depth\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Once we have a grid ID, we can then use the *pymt* *grid_* methods to get additional information\n",
+ "about each of the grids. Because this grid is uniform rectilinear (as returned by the\n",
+ "***grid_type*** method below), it is described by a set of methods that are only available\n",
+ "for grids of this type. These methods include:\n",
+ "* get_grid_shape\n",
+ "* get_grid_spacing\n",
+ "* get_grid_origin"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(\"Grid type: {0}\".format(cem.grid_type(2)))\n",
+ "print(\"Grid rank: {0}\".format(cem.grid_ndim(2)))\n",
+ "print(\"Grid shape: {0}\".format(cem.grid_shape(2)))\n",
+ "print(\"Grid spacing: {0}\".format(cem.grid_spacing(2)))\n",
+ "print(\"Grid origin: {0}\".format(cem.grid_origin(2)))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now that we know a little more about the grid, let's plot it with the current values of\n",
+ "water depth.\n",
+ "\n",
+ "Here I define a convenience function for plotting the water depth and making it look\n",
+ "pretty. You don't need to worry too much about it's internals for this tutorial.\n",
+ "It just saves typing later on."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_coast(spacing, z):\n",
+ " import matplotlib.pyplot as plt\n",
+ "\n",
+ " xmin, xmax = 0., z.shape[1] * spacing[0] * 1e-3\n",
+ " ymin, ymax = 0., z.shape[0] * spacing[1] * 1e-3\n",
+ "\n",
+ " plt.imshow(z, extent=[xmin, xmax, ymin, ymax], origin=\"lower\", cmap=\"ocean\")\n",
+ " plt.colorbar().ax.set_ylabel(\"Water Depth (m)\")\n",
+ " plt.xlabel(\"Along shore (km)\")\n",
+ " plt.ylabel(\"Cross shore (km)\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.empty(cem.grid_shape(2), dtype=float)\n",
+ "\n",
+ "cem.get_value(\"sea_water__depth\", out=z)\n",
+ "plot_coast(cem.grid_spacing(2), z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We have alredy set the incoming wave characteristics, but we've yet to add any sediment to\n",
+ "the system.\n",
+ "\n",
+ "From the list of input variables, can you tell which one might be the one we're looking\n",
+ "for?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for name in cem.input_var_names:\n",
+ " print(\"{0:70s} [{1}]\".format(name, cem.var_units(name)))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The one we want is *land_surface_water_sediment~bedload__mass_flow_rate*. Now have a look\n",
+ "at what sort of grid its defined on and how we could change its value."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.var_grid(\"land_surface_water_sediment~bedload__mass_flow_rate\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Notice that it's on the same grid as water depth. To add sediment, we need to:\n",
+ "1. allocate an array to hold sediment discharge values\n",
+ "2. set values of the sediment discharge array\n",
+ "3. pass this new sediment discharge array into CEM\n",
+ "\n",
+ "I've placed the sediment discharge in the horizontal center of the grid (column 100 of 200) and\n",
+ "along the bottom. The sediment will be routed in a straight line until it hits the coast.\n",
+ "\n",
+ "You don't need to do this, though. Feel free to add sediment in another location (or even multiple\n",
+ "locations!) or change the amount of sediment. Note that the CEM model is sensitive to the balance of\n",
+ "wave energy to sediment input. If you go too far from the defaults, you may get some \"interesting\"\n",
+ "results."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "qs = np.zeros_like(z)\n",
+ "qs[0, 100] = 750\n",
+ "cem.set_value(\"land_surface_water_sediment~bedload__mass_flow_rate\", qs)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now let's time step through the model. For every iteration of the for-loop we set the sediment\n",
+ "discharge (***set_value***), and then update the model to the next time (***update_until***)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for time in tqdm(range(3000)):\n",
+ " cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
+ " cem.update_until(time)\n",
+ "\n",
+ "cem.get_value('sea_water__depth', out=z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's have a look to see what sort of delta we've created after 3000.0 days."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cem.get_value('sea_water__depth', out=z)\n",
+ "plot_coast(cem.grid_spacing(2), z)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Exercise\n",
+ "\n",
+ "Now play with the CEM model on your own. To make things easier, I've placed all of the steps\n",
+ "to run CEM into a single cell below. Please feel free to modify!\n",
+ "\n",
+ "Some ideas\n",
+ "1. Modify wave energy vs sediment load\n",
+ "1. Change the incoming wave angle\n",
+ "1. Modify the river position so that it moves with time\n",
+ "1. Pick wave height and period from some probability density function\n",
+ "1. Add a second, or third, or fourth river\n",
+ "1. Increase the sediment load or have it vary with time\n",
+ "1. Make a movie of the evolving delta\n",
+ "\n",
+ "Anything else?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pymt.models\n",
+ "cem = pymt.models.Cem()\n",
+ "\n",
+ "args = cem.setup(number_of_rows=100, number_of_cols=200, grid_spacing=200.)\n",
+ "cem.initialize(*args)\n",
+ "\n",
+ "qs = np.zeros(cem.grid_shape(2), dtype=float)\n",
+ "qs[0, 100] = 750\n",
+ "\n",
+ "for time in tqdm(range(3000)):\n",
+ " cem.set_value(\"sea_surface_water_wave__height\", 1.5)\n",
+ " cem.set_value(\"sea_surface_water_wave__period\", 7.)\n",
+ " cem.set_value(\n",
+ " \"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\",\n",
+ " 0. * np.pi / 180.,\n",
+ " )\n",
+ "\n",
+ " cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
+ " cem.update_until(time)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.empty(cem.grid_shape(2), dtype=float)\n",
+ "cem.get_value('sea_water__depth', out=z)\n",
+ "\n",
+ "plot_coast(cem.grid_spacing(2), z)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def avulse_river(river_x, stddev=1.0, x_min=None, x_max=None):\n",
+ " river_x += np.random.normal(0.0, stddev)\n",
+ " if x_max is not None and river_x >= 200:\n",
+ " river_x = river_x - 200\n",
+ " if x_min is not None and river_x < 0:\n",
+ " river_x = 200 + river_x\n",
+ " return river_x"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "river_x = 100.0\n",
+ "river_i = []\n",
+ "for time in range(3000):\n",
+ " river_x = avulse_river(river_x, x_min=0.0, x_max=200.0)\n",
+ " river_i.append(int(river_x))\n",
+ "plt.plot(river_i)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Exercise\n",
+ "\n",
+ "### Couple models\n",
+ "\n",
+ "Instead of using constant wave characteristics, let's now couple the CEM component with a wave\n",
+ "gererator component.\n",
+ "\n",
+ "### Waves"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from pymt.models import Waves\n",
+ "\n",
+ "waves = Waves()\n",
+ "args = waves.setup(angle_asymmetry=0.2, angle_highness_factor=0.5)\n",
+ "waves.initialize(*args)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "angles = np.zeros(1000)\n",
+ "for day in range(1000):\n",
+ " waves.update()\n",
+ " angles[day] = waves.get_value(\"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "_ = plt.hist(angles * 180.0 / np.pi, bins=100)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now add the waves component to our coupling script."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pymt.models\n",
+ "\n",
+ "waves = pymt.models.Waves()\n",
+ "args = waves.setup(angle_asymmetry=0.2, angle_highness_factor=0.5)\n",
+ "waves.initialize(*args)\n",
+ "\n",
+ "cem = pymt.models.Cem()\n",
+ "args = cem.setup(number_of_rows=100, number_of_cols=200, grid_spacing=200.)\n",
+ "cem.initialize(*args)\n",
+ "\n",
+ "qs = np.zeros(cem.grid_shape(2), dtype=float)\n",
+ "qs[0, 100] = 500\n",
+ "\n",
+ "for time in tqdm(range(3000)):\n",
+ " waves.update()\n",
+ " angle = waves.get_value(\"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\")\n",
+ " \n",
+ " cem.set_value(\"sea_surface_water_wave__height\", 1.5)\n",
+ " cem.set_value(\"sea_surface_water_wave__period\", 7.0)\n",
+ " cem.set_value(\n",
+ " \"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\",\n",
+ " angle,\n",
+ " )\n",
+ "\n",
+ " cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
+ " cem.update_until(time)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = cem.get_value('sea_water__depth').reshape(cem.grid_shape(2))\n",
+ "plot_coast(cem.grid_spacing(2), z)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.8.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/lessons/pymt/cem_and_waves_and_hydrotrend.ipynb b/lessons/pymt/cem_and_waves_and_hydrotrend.ipynb
new file mode 100644
index 0000000..f38d15a
--- /dev/null
+++ b/lessons/pymt/cem_and_waves_and_hydrotrend.ipynb
@@ -0,0 +1,1193 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\u001b[33;01m➡ models: Avulsion, Plume, Sedflux3D, Subside, FrostNumber, Ku, Hydrotrend, Cem, Waves, ExponentialWeatherer, Flexure, FlowAccumulator, FlowDirectorD8, FlowDirectorDINF, FlowDirectorSteepest, LinearDiffuser, OverlandFlow, SoilMoisture, StreamPowerEroder, TransportLengthHillslopeDiffuser, Vegetation\u001b[39;49;00m\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "time = 0\n",
+ "time = 100\n",
+ "time = 200\n",
+ "time = 300\n",
+ "time = 400\n",
+ "time = 500\n",
+ "time = 600\n",
+ "time = 700\n",
+ "time = 800\n",
+ "time = 900\n",
+ "time = 1000\n",
+ "time = 1100\n",
+ "time = 1200\n",
+ "time = 1300\n",
+ "time = 1400\n"
+ ]
+ }
+ ],
+ "source": [
+ "import numpy as np\n",
+ "import pymt.models\n",
+ "\n",
+ "waves = pymt.models.Waves()\n",
+ "args = waves.setup()\n",
+ "waves.initialize(*args)\n",
+ "\n",
+ "cem = pymt.models.Cem()\n",
+ "args = cem.setup(number_of_rows=100, number_of_cols=200, grid_spacing=200.)\n",
+ "cem.initialize(*args)\n",
+ "\n",
+ "hydrotrend = pymt.models.Hydrotrend()\n",
+ "args = hydrotrend.setup(run_duration=100)\n",
+ "hydrotrend.initialize(*args)\n",
+ "\n",
+ "qs = np.zeros(cem.grid_shape(2), dtype=float)\n",
+ "qs[0, 100] = 750\n",
+ "qs_save = []\n",
+ "\n",
+ "for time in range(3000):\n",
+ " waves.set_value('sea_shoreline_wave~incoming~deepwater__ashton_et_al_approach_angle_asymmetry_parameter', .3)\n",
+ " waves.set_value('sea_shoreline_wave~incoming~deepwater__ashton_et_al_approach_angle_highness_parameter', .7)\n",
+ "\n",
+ " waves.update()\n",
+ " angle = waves.get_value(\"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\")\n",
+ " \n",
+ " cem.set_value(\"sea_surface_water_wave__height\", 2.0)\n",
+ " cem.set_value(\"sea_surface_water_wave__period\", 7.0)\n",
+ " cem.set_value(\n",
+ " \"sea_surface_water_wave__azimuth_angle_of_opposite_of_phase_velocity\",\n",
+ " angle,\n",
+ " )\n",
+ "\n",
+ " hydrotrend.update()\n",
+ " qs[0, 100] = hydrotrend.get_value(\"channel_exit_water_sediment~suspended__mass_flow_rate\")\n",
+ " \n",
+ " qs_save.append(qs[0, 100])\n",
+ " \n",
+ " cem.set_value('land_surface_water_sediment~bedload__mass_flow_rate', qs)\n",
+ " cem.update_until(time)\n",
+ " \n",
+ " if time % 100 == 0:\n",
+ " print(\"time = {0}\".format(time))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_coast(spacing, z):\n",
+ " import matplotlib.pyplot as plt\n",
+ "\n",
+ " xmin, xmax = 0., z.shape[1] * spacing[0] * 1e-3\n",
+ " ymin, ymax = 0., z.shape[0] * spacing[1] * 1e-3\n",
+ "\n",
+ " plt.imshow(z, extent=[xmin, xmax, ymin, ymax], origin=\"lower\", cmap=\"ocean\")\n",
+ " plt.colorbar().ax.set_ylabel(\"Water Depth (m)\")\n",
+ " plt.xlabel(\"Along shore (km)\")\n",
+ " plt.ylabel(\"Cross shore (km)\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAADnCAYAAAD7GCa6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfnUlEQVR4nO3de7gcVZ3u8e+bAIKACgYQBY0iiIhy9Rqd4XYcxlHBUZxBhoPKGbzgUccr+ni/oo+ijhckCBIVdXAUAVEUEWS8gYkgAUHxggrkECOiiNySvOePqoZOp/fu6t7d1dXN+3mefvbu2l2rfl3p/Pbaq1b9lmwTERHTad64A4iIiNFJko+ImGJJ8hERUyxJPiJiiiXJR0RMsQ3GHUBERNMceOCBXrVqVc/XLVu27Ju2D6whpIElyUdEdFi1ahVLly7t+TpJC2oIZ06S5CMiulgzJfcQJclHRHRYa7hjzdpxhzEUSfIREesxa9amJx8RMZUMSfIREdPKzph8RMRUS08+ImJKmfTkIyKmlu3MromImFa58BoRMc1y4TUiYnqlJx8RMdWcnnxExLRKWYOIiCmX4ZqIiCnlDNdEREyxzK6JiJhemV0TETHFUtYgImKKpaxBRMSUy3BNRMSUmqbhmnnjDiAiomnsoiff69GLpO0lnS/pSklXSHpFuX1LSedKurr8usWo3kuSfEREF8NI8sBq4NW2Hwk8ATha0i7AMcB5tncEziufj0SGayIiOgzrwqvtFcCK8vubJV0JPAg4CNinfNkS4ALg9XM+YBdJ8hERHfoYk18gaWnb88W2F3d7oaSFwB7ARcA25S8AbK+QtPWcAp5FknxERIc+kvwq23v3epGkzYAvA6+0/RdJc4ywuiT5iIhOHt4USkkbUiT4U21/pdx8g6Rty178tsDKoRysi1x4jYjo0CpQ1uvRi4ou+0nAlbaPa/vRmcAR5fdHAGcM/U2U0pOPiOgwxNo1i4DDgeWSLi23vRE4FjhN0pHA74BDhnGwbpLkIyI6eEiLhtj+HjDTAPz+cz5ABUnyERFdpKxBRMSUyqIhERFTzEOcXTNuSfIREV2kJx8RMaWmqQplknxERIcsGhIRMcWmaY3Xkd3x2oQ6yhERgxrGHa9NMMqyBmOvoxwRMYhhLRrSBCNL8rZX2P5J+f3NQHsd5SXly5YAB48qhoiIQU1LT76WMflB6ihLOgo4CmDTTTfda+edd64j1IiYcMuWLVtle6u5tLE2F16rG7SOcll4fzHAXnvv7R9cdPHogoyIqXGvDeb/dhjtTMpwTC8jTfLjrqMcETGIzJOvoEId5WOpWEfZNrevno4/nSJiAqSsQSVDq6O8dkhlPyMiqkiBsgqGWUfZwB1rk+Qjoh7TdDNU7niNiOgwrEVDmmAikrxt7pySEx4RkyHDNTWyyYXXiKiNmZw7WnuZiCS/1ubW1WvGHUZE3ENk0ZCIiCmX4Zoamem5CBIRzZfZNTWbpjoSEdF8tqdm2vZEJHlDZtdERG3Sk4+ImGbOmHytpmm9xYhovvTkazZNd59FRPOlCmXNjLl9TebJR0Q9bHPHlNyAORFJPiKibunJ12it4dYp+a0aEc2XO15rZmcKZUTUJ2PyNVtLZtdERJ2yaEhExNTKFMqapdRwRNRpmqZtT0iSz3BNRNQnPfmaZY3XiKiVp2dMft64A4iIaJpWT77XowpJJ0taKenytm1vk3SdpEvLx9NG9V569uQlPRH4N+ApwLbArcDlwNnA52z/eVTBtay1ue3O1aM+TETEXYbYkz8F+BjwmY7tH7L9gV47S5oH7AY8kCL/XmH7hqoHnzXJS/oGcD1wBvBuYCWwMbATsC9whqTjbJ9Z9YCDKC6CTMefThHRfMO8Gcr2hZIW9rufpB2A1wMHAFcDf6DMv5L+BpwALLE961h2r5784bZXdWz7K/CT8vFBSQv6Db5fWRkqIupkalk05GWS/jewFHi17T91/PxdwPHAi+x1/6yQtDXwPOBwYMlsB5l1TL4zwUu6j6QtW49ur2l77VjHoSIiBtXHmPwCSUvbHkdVPMTxwA7A7sAK4IPrxWAfavvCzgRf/myl7Q/bnjXBQ8XZNZJeBLyDYjyodUADD5tlt1OYwzhUuyz/FxG1qr5oyCrbe/fdfNuYuqQTga/N9FpJ84F/AhbSlrNtH1flWFWnUL4GeNRMvfZuBh2H6toWGa6JiJqNcJ68pG1tryifPotiMstMzgJuA5YDfSfCqkn+V8Df+m18Br3GoSIixm9Is2skfQHYh2Jo51rgrcA+knan6MNeA7xolia2s/2YQY9fNcm/AfiBpIuA21sbbb+8z+MdD7yT4o29k2Ic6oXdXliObR0FsOlWD0hPPiLqU0yvGVJTPrTL5pP6aOIbkp5q+1uDHL9qkj8B+A4D/rnQ0s84lO3FwGKA++/wSCfJR0StmlPW4EfA6eV8+TsBAbZ9nyo7V03yq22/asAA79LnONRd1mJuXZ3l/yKiJqZJSf6DwBOB5d1m2vRSNcmfXw6fnMW6wzU3zrTDEMahIiLGpzm1a64GLh8kwUP1JP+88usb2rbNOoVyCONQd1mbUsMRUSs3qSe/ArigrEDQ3ske6hTKh3W542rjyiHOkZ0qlBFRs+b05H9TPjYqH32pmuRPom0WjKRNgTOB/fs94CCMs8ZrRNSnuOV13FEAYPvtc9m/aqnh6yQdDyBpC+Bc4HNzOXBERHO56Mn3eoyQpMWSHj3DzzaV9EJJh/Vqp1JP3vabJb1P0ieBvYBjbX+5v5AHN01LcUXEBGjG7JpPAG8uE/3l3F2FckfgPsDJwKm9GulVavif255eDLy5/GpJ/2z7K4PF3h/bufAaEfUa85i87UuB50raDNibu9fzuNL2z6u206sn/4yO55cAG5bbDdSS5NeSnnxE1Gz8PXkAbP8VuGDQ/WdN8rZfMGjDERETa4hlDcat13DNm4CPz1RETNJ+wL1tz1ieYBjs3PEaETVrSE9+rnoN1ywHvibpNoqVoNoH/ncHvg28Z6QRUpzrO9YkyUdEjZozT35Oeg3XnEGxjuuOwCKKgf+/UEyfPMr2raMPsVVPfjpOeERMANOYJC9pJ+C1wENYd9GQ/arsX3UK5dUU9RMiIu4ZmjNc8yXgk8CJQN9DGlXveB0rZ/m/iKhVo2rXrLZ9/KA7T0iSzxTKiKhRA8oaSNqy/PYsSS8FTqdiFeB2k5HkSU8+Imo2/jH5ZRS/blQ+f23bz2atAtyuUpIvB/6PB7axvaukxwDPtP2u6vFGREyQ8d/x+lAoKv7avq39Z/1UAa7akz+R4rfICeXBL5P0eaCWJG/DnZknHxF1caPG5H8A7FlhW1dVk/y9bV8sqX3b6or7zt0U3X0WERNizD15SQ8AHgRsImkP7h62uQ9w76rtVE3yqyTtQDEOhKTnUKxWUp8k+Yio0/h78v8APB/YDmhfBepm4I1VG6ma5I8GFgM7S7qOYpWSnnWMIyImUgNm19heAiyR9Oy5lHbvmeQlzQP2tn1AuSLUPNs3D3rAgRhIqeGIqM3oFwXpwwWS/hN4MkU2/B7wDtt/rLJzzyRve62klwGn2b5lTqEOKmPyEVGnZiwa0vJF4ELg2eXzw4D/Ag6osnPV4ZpzJb2mbPiuRF91Mv5QZCHviKhTc3ryW9p+Z9vzd0k6uOrOVZN8axHvo9u2VZ6MHxExcZrTkz9f0r8Cp5XPnwOcXXXnqgXKHjpAYMOT4ZqIqFOzcs6LgFdRVP81MB+4RdKrANu+z2w7V73jdUPgJcDflZsuAE6wfeeAQffHwJ25GSoiatSQnrztzeeyf9XhmuMp1nb9RPn88HLb/5lpB0knA08HVtretdy2JcW4/kLgGuC5M606ta5G/VaNiGnXrHryorjY+lDb75S0PbCt7Yur7D+v4nEea/sI298pHy8AHttjn1OAAzu2HQOcZ3tH4LzyeURE86x170c9PgE8EXhe+fyvwMer7ly1J79G0g62fwUg6WH0KF5v+0JJCzs2HwTsU36/hGLY5/U9j96AGxMi4h6mIT154PG295R0CYDtP0naqOrOVZP8aymu8P6aon7CQ4AX9B1qUcVyBYDtFZK2numFko4CjgJg0y0ga7xGRG0adTPUnZLmc3dZma2Ayr3eqrNrzivXeX0ERZK/yvbtPXabE9uLKUopoPs/2GSN14ioS7NGD/6TYsGQrSW9m2IK5Zuq7tzPoiF7UVww3QDYTRK2P9PH/gA3SNq27MVvC6zsc/+IiHo0Z3bNqZKWAftTdLIPtn1l1f2rTqH8LLADcCl3j8Ub6DfJnwkcARxbfj2j0l7N+q0aEdOuQWUNJD0a2JmiU3xlPwkeqvfk9wZ2sasPUkn6AsVF1gWSrgXeSpHcT5N0JPA74JBqrRmyaEhE1Gb8Y/KS7kvREd4euIyiF/9oSb8DDrL9lyrtVE3ylwMPoI8a8rYPneFH+1dt4+7GSE8+Iuo1pJ78HO4ZeiewFNjP9tpyv/nAe4F3A/+3yvFnnScv6SxJZwILgJ9J+qakM1uPam8xImLCtG6G6vWo5hQGu2foAOCYVoIHsL2GYsGQShUooXdP/gNVGxqpZtWRiIh7giHlnDncM3SH7fWWWbW9WlLl2Y2zJnnb3wUoFwu5tawtvxPFRYBvVD3IUCTJR0Rdqi/kvUDS0rbni8vp371UuWdo4461XVsE3KtKcFB9TP5C4CmStqD402Ip8C/UtQRgevIRUbdqwzGrbO89oghWsO7aru3+X9VGqiZ52f5bOSvmo7bfL+nSqgeJiJg4o51d0/OeIdv7DuNAlZO8pCdS9NyPLLfNH0YAlWSN14io22jnyQ92z9AAqib5VwBvAE63fUVZoOz8UQXVVYZrIqIuQxwiHu49Q/2rWrvmQopx+dbzXwMvH1VQXQLIzVARUa8hDdfM5Z6hspb8drZ/P+jxq9aTj4i452iVNRhzPfmyysBX59JGPwXKxid3vN4z3HrH3d9vslHxfJON7v7ZJpVLaDfPNL+3OrTOX53nqTmlhn8k6bG2fzzIzpOR5PHd/8gbzYc7OoZuNhriNeCZ2m5tbz/+bD8bRWwzqRLzMNue6f12Ozf9tt/S+vduT47t37eO1es4nbHOdG6qtDWozvfY6721x1WX2c5NXcdvHXOmzwR0P0+j+rdrSIEyYF/gxZKuAW6hmCdv24+psnPVKpTvB94F3AqcA+wGvNL25waJuG/tJ7vbB6B92/x56/b65/cYkVqzdv19Zmq78/lsP2vf1oqh27GqxNv6+WxxVo25Xbf22o/f7VhV33/r+/b3XqX9qlrtt9YZ6Gy/s+3W6/v592w/PzN9jvr9t6lqts90Zzy9PsPdzkm31890blqvnenczvbzXuet02wJfrbX97tfL9VvhqrDP85l56pj8k8tK549HbgW2IlitaiIiOm0dm3vRw1s/5aiEuV+5fd/o4/rqVWHazYsvz4N+ILtG4uLvg3U2UOo0rsa9Xh/e/u94pstlmHH2a29UR9jktqf7d+tVxzDVOXfqd/PTT/xtl470z6z/XySr6U1pCcv6a0U5d4fAXyaIh9/DlhUZf+qSf4sSVdRDNe8tFxj8Lb+w42ImACtKpTN8CxgD+AnALavl7R51Z2rzpM/RtL7gL/YXiPpFooqahHD88z3r7/tzNfVH8eotd7nNL63qdGoMfk7bFtSayHvTfvZudK4jqRDgNVlgn8TxZ8KD+w71IiISTG8evJzdZqkE4D7Sfp34NvAp6ruXHW45s22vyTpycA/UNSZPx54fL/RDkV7j68JvaFuPVCoL7ZRH7+O9zfTMYb5b915jDo/O93e3zjjmU0T/tKY6fPQMurYGrTGq+0PSPpfwF8oxuXfYvvcqvuryrKtki6xvYek9wLLbX++tW3gyPugez/AvO6G2V80l3/02T5Qvdrt9WGs0sZcjPr4VdqfyzGqtj/o8UYd/zBjaBlHcu0nxlHEN+jnANaP55L3L5tr+V9tuJXZosKI9B9OmvOxesYivc/263ttm3H/ikn+a8B1FEtO7UVxAfZi27v1H3L/9ECZo/rYoZ8P4VySwCAfzPZ25joGPdfjj7v9QY9R9XjDTByDGOV76+c4w+iozKX9URxzthiGkeQ3WGDuVyHJ//HkOpL8T2zv2bHtsqHeDAU8l2KNwg/Yvqmsf5x58hExvcY8u0bSS4CXAg+TdFnbjzYHvl+5nSo9+fKAuwFPKZ/+j+2fVj3IXPXdk283Uy9j2D2JOrTey7Bir+PcjKqn3etYdf371vX56jxOkz+/o/7LoZe3M5ye/ObP6P3Cm04ZWU9e0n2BLYD3su5C3zfbvrFyOxWHa14B/DvwlXLTsyjWMvxo5YjnYE5JPiLuWYaV5Dd7eu8X/nnJyIdrWsp1YDduPbf9uyr7VR2uORJ4vO1byoO9D/ghUEuSj4ioVYPWlZb0DIq1Xh9IsUzgQ4ArgUdV2b/y8n9AewWgNay/gnhlZTW1m8t2Vtf1mzAiorKGTKGkKA75BODb5SzHfYGZFiJZT9Uk/2ngIkmnl88PBk7qK8z17Wt71RzbiIgYvmaVNbjT9h8lzZM0z/b55WhKJVXLGhwn6QLgyRQ9+BfYvmSweCMiJkBzevI3SdqMYgnWUyWtBFZX3blnkpc0D7jM9q6UBXKGwMC3yloMJ9hePKR2IyKGozk9+YMoCkL+B3AYcF/gHVV37pnkba+V9FNJD656NbeCRWUlta2BcyVdVS4WfhdJR0E5p+a+QzpqREQVDVg0RNIrKebDX2K7dU10Sb/tVB2T3xa4QtLFFMtPAWD7mf0esNzv+vLrynKc/3EUf4q0v2YxsBjKKZQREXUa/+ya7YCPADuXN0P9gCLp/7CfefJVk/zb+4+vu7JM5jzbN5ffP5U+/vSIiKjFmIdrbL8GQNJGFIuGPAl4IXCipJts71KlnVmTvKSHA9vY/m7H9r+jqGUziG2A08uVpTYAPm/7nAHbiogYjSoDCPX8HtgEuA/FwPV9geuB5VV37tWT/zDwxi7b/1b+rMJ9v+uy/WuKhcAjIhrKMK9CBh/y+uHtJC2muOHpZuAiiuGa42z/qZ92eiX5hbYv69xoe6mkhf0cKCJiYoixJ3ngwcC9gKspRk6uBW7qt5FeSX7jWX62Sb8Hi4iYGPMrXHi9c3SHt32ginHtR1GMx78a2FXSjRQXX99apZ1ey//9uFxuah2SjgSW9RlzRMSEcDEm3+sx6igKlwNfB75BMbtmB+AVVdvo1ZN/JcVF0sO4O6nvDWxEUYkyImL6iFqS+KwhSC+n6MEvovib4fsUhSFPZlgXXm3fADypLIiza7n5bNvfGSToiIiJUWVMfrQWAv8N/IftFYM2UrV2zfnA+YMeJCJi4ow5ydt+1TDaqXozVETEPUcDhmuGJUk+ImI9rja7ZgIkyUdEdKo6T34CJMlHRHST4ZqIiCmWJB8RMa0q1q6ZAEnyERGdhji7RtI1FEXG1gCrbe89lIYrSpKPiOhmg6HOrtnX9qphNlhVknxERKcpmiffq0BZRMQ9UDkm3+sBCyQtbXsc1b0xviVp2Qw/H6n05CMiuqnWk19VYYx9ke3rJW0NnCvpKtsX9thnaNKTj4jo1LoZqndPvifb15dfVwKnA48bXeDrS5KPiOhm/trejx4kbSpp89b3wFOBy0cc+ToyXBMRsZ6hLQqyDcWaHFDk28/bPmcYDVeVJB8R0WlItWts/xrYbc4NzUGSfEREN1MyhTJJPiKim5Q1iIiYUlN0M1SSfETEeqZn0ZCxTKGUdKCkn0v6paRjxhFDRMSs5N6PCVB7kpc0H/g48I/ALsChknapO46IiBkN8WaocRtHT/5xwC9t/9r2HcAXgYPGEEdExMymJMmPY0z+QcDv255fCzy+80VlIZ9WMZ/beXu9d4l1sQAYS6nQDk2IowkxQDPiaEIM0Iw4mhADwCPm3sTkDMf0Mo4kry7b1jubthcDiwEkLa270H6nJsTQlDiaEENT4mhCDE2JowkxtOKYeyNMTE+9l3Ek+WuB7duebwdcP4Y4IiJmltk1A/sxsKOkh0raCPhX4MwxxBER0V1rnvwUzK6pvSdve7WklwHfBOYDJ9u+osdui0cfWU9NiAGaEUcTYoBmxNGEGKAZcTQhBhhKHJNzYbUX2dPxRiIihkWbb2L22qH3C797xbImXIeYTe54jYjolAuvERFTLhdeR68p5Q8kXSNpuaRLhzI9q/pxT5a0UtLlbdu2lHSupKvLr1uMIYa3SbquPB+XSnraiGPYXtL5kq6UdIWkV5Tb6z4XM8VR2/mQtLGkiyX9tIzh7eX2us/FTHHU+tkojzlf0iWSvlY+H8K5qHDRdUIuvDY2yTew/MG+tnevefztFODAjm3HAOfZ3hE4r3xedwwAHyrPx+62vz7iGFYDr7b9SOAJwNHlZ6HuczFTHFDf+bgd2M/2bsDuwIGSnkD952KmOKDezwbAK4Ar257P/VykrEEt7vHlD8oV3W/s2HwQsKT8fglw8BhiqJXtFbZ/Un5/M8V/6AdR/7mYKY7auPDX8umG5cPUfy5miqNWkrYD/gn4VNvm4ZyL9ORHrlv5g1r/Q7Ux8C1Jy8pyC+O0je0VUCQdYOsxxfEySZeVwzkjHRpoJ2khsAdwEWM8Fx1xQI3noxyeuBRYCZxreyznYoY4oN7PxoeB1wHtA+jDORfpyY9cpfIHNVlke0+KoaOjJf3dmOJoiuOBHSj+TF8BfLCOg0raDPgy8Erbf6njmBXjqPV82F5je3eKu8UfJ2nXUR6vzzhqOxeSng6stL1s+I2TnnwNGlP+wPb15deVwOkUQ0njcoOkbQHKryvrDsD2DeV/8LXAidRwPiRtSJFYT7X9lXJz7eeiWxzjOB/lcW8CLqC4ZjK2z0V7HDWfi0XAMyVdQzGcu5+kzzGUc1EuGtLrMQGanOQbUf5A0qaSNm99DzwVxloR80zgiPL7I4Az6g6g9R+o9CxGfD4kCTgJuNL2cW0/qvVczBRHnedD0laS7ld+vwlwAHAV9Z+LrnHUeS5sv8H2drYXUuSH79j+N4ZxLqbowmtj58kPWP5gFLYBTi/+f7MB8Hnb59RxYElfAPYBFki6FngrcCxwmqQjgd8Bh4whhn0k7U4xfHYN8KJRxkDRYzscWF6OAQO8kZrPxSxxHFrj+dgWWFLOPpsHnGb7a5J+SL3nYqY4PlvzZ6Ob4XwuJmQ4ppeUNYiI6KAt7mX2f0DvF375dylrEBExcSbowmovSfIREd1MyJh7L0nyERGdxMTMnuklST4iopsM10RETKvJmSLZS5PnyccYSXqWJEvauW3bQrVVo6wxln1aFQZrONbBkt5Sfn+KpOcM2M5WkmqZahsjkjteY8odCnyP4iaTiVbO5a7qdcAn5npM238AVkhaNNe2Ygym6GaoJPlYT1mbZRFwJDMk+bKe+KdV1Nm/RNK+5fbnS/qKpHPKet7vb9vnSEm/kHSBpBMlfaxLu3/fVov8ktbdxsBmkv5b0lWSTi3vPkXS/uXrlpcFse5Vbr9G0lskfQ84RNIOZUzLJP1P+18obcfeCbjd9qouP3tn2bOfV7b9Hkk/lLRU0p6SvinpV5Je3LbbV4HDqp31aJwpKWuQMfno5mDgHNu/kHSjpD1bJXbbHA1g+9FlwvxWmSShKE61B0XN8Z9L+iiwBngzsCdwM/Ad4Kddjv0a4Gjb3y9/2dxWbt8DeBRF/aLvA4tULOByCrB/GetngJdQVCYEuM32kwEknQe82PbVkh5P0Vvfr+PYi4DO90n5i+q+wAtsu/z98nvbT5T0oTKGRcDGwBXAJ8tdlwLv6vIeo+kmaDiml/Tko5tDKQo+UX49tMtrngx8FsD2VcBvgVaSP8/2n23fBvwMeAhFoarv2r7R9p3Al2Y49veB4yS9HLif7dXl9ottX1sWvroUWAg8AviN7V+Ur1kCtFcI/S+46y+TJwFfKssRnEBxW36nbYE/dGx7cxnHi7zu7eGtOkrLgYts31wO0dzWqulCURjrgTO8z2i6KRmuSU8+1iHp/hQ93F0lmaJukCW9rvOlszRze9v3ayg+Z7O9/i62j5V0NvA04EeSDphDm7eUX+cBN5VlcWdzK0WPvd2Pgb0kbWm7ffGUVjxrO2Jby93/rzYu24xJlJ58TKnnAJ+x/RDbC21vD/yGoufe7kLK8eZymObBwM9nafdi4O8lbSFpA+DZ3V4kaQfby22/j2K4Y72x8zZXAQslPbx8fjjw3c4XlTXffyPpkPIYkrRbl/auBB7ese0cioJXZ7ddH6hqJ8ZbsTQGlQuvMcUOpaiZ3+7LwPM6tn0CmC9pOcWwyPNt384MbF8HvIdiJaVvUwzj/LnLS18p6XJJP6XoBX9jljZvA15AMQyznKIX/ckZXn4YcGTZ7hV0X0ryQmCP1kXdtuN8iaI2+pkqyupWtS9wdh+vjyaZkimUqUIZtZG0me2/lj350ynKR3f+QhkrSR8BzrL97SG0dSFwkO0/zT2yqJO2nm/+ZdPeL/zYzT2rUEo6EPgIxdDnp2wfO5QgK0pPPur0tvLC5+UUQ0BfHXM83bwHuPdcG5G0FXBcEvyEGtJwTXmPxscplg7dhWLtgV1GG/y6cuE1amP7NeOOoRfbNzCEFcjKmTZN/CUWVQ1nOOZxwC9t/xpA0hcphgp/NozGq0iSj4joptqF1QXl/Roti20vbnv+IOD3bc+vBR4/hOgqS5KPiOhU/cLqqh5j8t2m+dZ6ITRJPiKim+FMkbwW2L7t+XYUd23XJkk+IqLT8BYN+TGwo6SHAtdR1ILqnI48UknyERHdDOHCq+3Vkl4GfJNiCuXJtq+Yc8N9SJKPiFjP8O5otf114OtDaWwASfIREd1MyB2tvSTJR0R0at0MNQWS5CMiupmQRUF6Se2aiIgO5fq8Cyq8dJXtA0cdz1wkyUdETLEUKIuImGJJ8hERUyxJPiJiiiXJR0RMsST5iIgp9v8BbszGTH0REiMAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "