{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A notebook to generate a river forcing. Constant Fraser river outflow 2000m^3/s.\n",
    "\n",
    "Based on /data/nsoontie/MEOPAR/NEMO-forcing/rivers/Fraser_only_cnst.nc\n",
    "\n",
    "Uses a runoff temperaure of 14 deg to avoid runaway temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import netCDF4 as nc\n",
    "import numpy as np\n",
    "\n",
    "from salishsea_tools import nc_tools\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3D constant forcing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First look at the constant Fraser river forcing file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "river = nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/rivers/rivers_Fraser_only_cnst.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KeysView(OrderedDict([('nav_lat', <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 nav_lat(y, x)\n",
      "unlimited dimensions: \n",
      "current shape = (898, 398)\n",
      "filling off\n",
      "), ('nav_lon', <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 nav_lon(y, x)\n",
      "unlimited dimensions: \n",
      "current shape = (898, 398)\n",
      "filling off\n",
      "), ('time_counter', <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 time_counter(time_counter)\n",
      "    units: non-dim\n",
      "unlimited dimensions: time_counter\n",
      "current shape = (1,)\n",
      "filling off\n",
      "), ('rorunoff', <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 rorunoff(time_counter, y, x)\n",
      "    _Fillvalue: 0.0\n",
      "    _missing_value: 0.0\n",
      "    _units: kg m-2 s-1\n",
      "unlimited dimensions: time_counter\n",
      "current shape = (1, 898, 398)\n",
      "filling off\n",
      "), ('rodepth', <class 'netCDF4._netCDF4.Variable'>\n",
      "float32 rodepth(y, x)\n",
      "    _Fillvalue: -1.0\n",
      "    missing_value: -1.0\n",
      "    units: m\n",
      "unlimited dimensions: \n",
      "current shape = (898, 398)\n",
      "filling off\n",
      ")]))\n"
     ]
    }
   ],
   "source": [
    "nc_tools.show_variables(river)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "run_off_3D = river.variables['rorunoff'][0,:,:]\n",
    "depth_3D =river.variables['rodepth'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inds = np.where(run_off_3D != 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.31223202,  4.62446404,  2.31223202], dtype=float32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_off_3D[inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.24893\n"
     ]
    }
   ],
   "source": [
    "total_run_off_3D = np.sum(run_off_3D[inds])\n",
    "print(total_run_off_3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(depth_3D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dishcharge distrubuted across three different points.\n",
    "\n",
    "Plan: add these up and distrubute across right end of my domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Towards 2D files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load bathymetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fB=nc.Dataset('/data/nsoontie/MEOPAR/2Ddomain/grid/bathy2D_36.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iniitialize run off arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d = fB.variables['Bathymetry'][:]\n",
    "X = fB.variables['x'][:]\n",
    "Y=fB.variables['y'][:]\n",
    "ymax, xmax = d.shape\n",
    "runoff = np.zeros((ymax, xmax))\n",
    "run_depth = -np.ones((ymax, xmax)) \n",
    "run_temp = -999*np.ones((ymax, xmax)) # -999 means missing data and will use surface temp\n",
    "\n",
    "TEMP = 14 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "run_off_depth = 3 # fresh water will be added over three metres"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add river to second to last grid cell. Divide discgarge in spanwise, but exlude first and last."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.25\n"
     ]
    }
   ],
   "source": [
    "#depth at last grid cell?\n",
    "print(d[2,-2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "js=1; je=ymax-2\n",
    "\n",
    "runoff[js:je+2,xmax-2] = total_run_off_3D/(je-js)\n",
    "run_depth[js:je+2,xmax-2]= run_off_depth\n",
    "run_temp[js:je+2,xmax-2]= TEMP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f07249690b8>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAET5JREFUeJzt3X+s3XV9x/HnyxbMmD/QkKACpji7BRJRdBY2f3Dc3LyS\naRf/GHbqJsgkS9ptWZyIS+Qas0W2kDlTxeqAODLpEiWuJPyQTU9GnENxgKhU2zC2Vib+wh/zR2zD\ne3+c03J2bO/5Xu739p5zvs9HctPvj8/9nvf79vZ1Pvdzzvc2VYUkab49bq0LkCStPsNekjrAsJek\nDjDsJakDDHtJ6gDDXpI6YGLYJ7kmyUNJ7p0w7oVJDiZ5TXvlSZLa0GRmfy2wsNSAJOuAK4BbgLRQ\nlySpRRPDvqpuBx6eMGwb8FHgm20UJUlq14rX7JOcAmwGrhoe8pZcSZoybbxA+x7gbTX4vQvBZRxJ\nmjrrW7jGC4CdSQBOAl6Z5EBV7RodlMQZvyQ9BlW14kn0imf2VfWsqjq9qk5nsG7/h+NBPzJ2bj8u\nv/zyNa/B/uzN/ubvoy0TZ/ZJrgfOA05Ksg+4HDhuGN47WqtEkrRqJoZ9VW1perGqunBl5UiSVoN3\n0Lak1+utdQmrap77m+fewP40kDbXhJZ8oKSO1WNJ0rxIQk3DC7SSpOln2EtSBxj2ktQBhr0kdYBh\nL0kdYNhLUgMHD8K9S/6vHtPNsJekBq67Ds46a62reOwMe0lq4Ec/WusKVsawl6QOMOwlqQMMe0nq\nAMNekhqY9V/tZdhLUgcY9pLUQGb8f9c27CWpAZdxJElTz7CXpAZcxpGkDnAZR5I09Qx7SWrAZRxJ\n6oC5X8ZJck2Sh5Ic8Tc5J3ldknuSfCHJp5PM8C8BlaT51GRmfy2wsMT5+4GXVtVZwLuAD7ZRmCRN\nk7lfxqmq24GHlzj/mar63nD3DuDUlmqTJLWk7TX7NwE3tXxNSdIKrW/rQkleBlwEvOhoYxYXFw9v\n93o9er1eWw8vSXOh3+/T7/dbv26qwUvMSTYAN1bVc45y/izgBmChqvYeZUw1eSxJmkbvex9s3Xrs\n35WThKpa8SsGK17GSfJMBkH/+qMFvSRpbU1cxklyPXAecFKSfcDlwHEAVbUDeAfwFOCqDF6uPlBV\nm1atYknSsk0M+6raMuH8xcDFrVUkSWqdd9BKUgcY9pLUAYa9JHWAYS9JHWDYS1IHGPaS1AGGvSR1\ngGEvSR1g2EtSBxj2ktQBhr0kdYBhL0kdYNhLUgcY9pLUAYa9JHWAYS9JHWDYS1IHGPaS1AGGvSR1\ngGEvSQ1UrXUFK2PYS1IHGPaS1ECy1hWsjGEvSQ3M/TJOkmuSPJTk3iXGvDfJniT3JDm73RIlSSvV\nZGZ/LbBwtJNJzgeeXVUbgTcDV7VUmyRNjblfxqmq24GHlxjyauDDw7F3ACcmObmd8iRpOsz9Mk4D\npwD7Rvb3A6e2cF1JUkvWt3Sd8R9wjvgcuLi4eHi71+vR6/VaenhJWl3Hahmn3+/T7/dbv26qwc8m\nSTYAN1bVc45w7gNAv6p2Dvd3A+dV1UNj46rJY0nSNNq+HbZtO/bLOUmoqhU/1bSxjLML+D2AJOcC\n3x0PeknS2pq4jJPkeuA84KQk+4DLgeMAqmpHVd2U5Pwke4EfAheuZsGStBZm/d04E8O+qrY0GLO1\nnXIkaTrN+iq0d9BKUgcY9pLUwKwv4xj2ktSAyziSpKln2EtSAy7jSFIHuIwjSZp6hr0kNeAyjiR1\ngMs4kqSpZ9hLUgMu40hSB7iMI0maeoa9JDXgMo4kdYDLOJKkqWfYS1IDLuNIUge4jCNJmnqGvSQ1\n4DKOJHWAyziSpKln2EtSB0wM+yQLSXYn2ZPk0iOcPynJLUnuTvLFJG9clUolSY/ZkmGfZB2wHVgA\nzgS2JDljbNhW4K6qeh7QA65Msn4VapUkPUaTZvabgL1V9UBVHQB2ApvHxvwP8KTh9pOAb1fVwXbL\nlCStxKQZ+CnAvpH9/cA5Y2M+BHwyyYPAE4Hfaa88SVIbJoV9kzcbvR24u6p6SX4BuC3Jc6vqB+MD\nFxcXD2/3ej16vd4ySpWk+dfv9+n3+61fN7XEm0eTnAssVtXCcP8y4JGqumJkzE3AX1TVp4f7/wJc\nWlV3jl2rlnosSZpm27fDtm3H/v32SaiqFd/SNWnN/k5gY5INSY4HLgB2jY3ZDbx8WNTJwC8B96+0\nMElSe5Zcxqmqg0m2ArcC64Crq+q+JJcMz+8A/hK4Nsk9DJ483lpV31nluiVJyzDxLZJVdTNw89ix\nHSPb3wJe1X5pkqS2eAetJHWAYS9JHWDYS1IHGPaS1AGGvSR1gGEvSR1g2EtSBxj2ktQBhr0kdYBh\nL0kdYNhLUgcY9pLUAYa9JHWAYS9JHWDYS1IHGPaS1AGGvSR1gGEvSR1g2EtSBxj2ktQBhr0kdYBh\nL0kdYNhLUgdMDPskC0l2J9mT5NKjjOkluSvJF5P0W69SkrQi65c6mWQdsB14OfA14HNJdlXVfSNj\nTgTeB7yiqvYnOWk1C5YkLd+kmf0mYG9VPVBVB4CdwOaxMb8LfKyq9gNU1bfaL1OStBKTwv4UYN/I\n/v7hsVEbgacm+VSSO5O8oc0CJUkrt+QyDlANrnEc8Hzg14ETgM8k+feq2jM+cHFx8fB2r9ej1+s1\nLlSSuqDf79Pv91u/7qSw/xpw2sj+aQxm96P2Ad+qqh8DP07yr8BzgSXDXpL0s8Ynwu985ztbue6k\nZZw7gY1JNiQ5HrgA2DU25p+AFydZl+QE4Bzgy61UJ0lqxZIz+6o6mGQrcCuwDri6qu5Lcsnw/I6q\n2p3kFuALwCPAh6rKsJekKZKqJsvyLTxQUsfqsSSpbdu3w7ZtcKxjLAlVlZVexztoJakDDHtJ6gDD\nXpI6wLCXpA4w7CWpAwx7SeoAw16SOsCwl6QOMOwlqQMMe0nqAMNekjrAsJekDjDsJakDDHtJ6gDD\nXpI6wLCXpA4w7CWpAwx7SeoAw16SOsCwl6QOMOwlqQMMe0nqAMNekjpgYtgnWUiyO8meJJcuMe6F\nSQ4meU27JUqSVmrJsE+yDtgOLABnAluSnHGUcVcAtwBZhTolSSswaWa/CdhbVQ9U1QFgJ7D5COO2\nAR8FvtlyfZKkFkwK+1OAfSP7+4fHDktyCoMngKuGh6q16iRJrVg/4XyT4H4P8LaqqiRhiWWcxcXF\nw9u9Xo9er9fg8pLUHf1+n36/3/p1U3X0PE9yLrBYVQvD/cuAR6rqipEx9/NowJ8E/Aj4g6raNXat\nWuqxJGmabd8O27bBsY6xJFTVil8LnTSzvxPYmGQD8CBwAbBldEBVPWukqGuBG8eDXpK0tpYM+6o6\nmGQrcCuwDri6qu5Lcsnw/I5jUKMkaYUmzeypqpuBm8eOHTHkq+rCluqSJLXIO2glqQMMe0nqAMNe\nkjrAsJekDjDsJakDDHtJ6gDDXpI6wLCXpA4w7CWpgUO/E2dWf8WXYS9JHWDYS1IHGPaS1IDLOJKk\nqWfYS1IHGPaS1IDLOJKkqWfYS1IHGPaS1IDLOJKkqWfYS1IHGPaS1IDLOJKkqWfYS1IHNAr7JAtJ\ndifZk+TSI5x/XZJ7knwhyaeTnNV+qZK0duZ+GSfJOmA7sACcCWxJcsbYsPuBl1bVWcC7gA+2Xagk\n6bFrMrPfBOytqgeq6gCwE9g8OqCqPlNV3xvu3gGc2m6ZkqSVaBL2pwD7Rvb3D48dzZuAm1ZSlCRN\nm1lfxlnfYEzj1pK8DLgIeNGRzi8uLh7e7vV69Hq9ppeWpE7o9/v0+/3Wr5ua8DSV5FxgsaoWhvuX\nAY9U1RVj484CbgAWqmrvEa5Tkx5LkqbVlVfCW94CP/kJPP7xx+5xk1BVWel1mizj3AlsTLIhyfHA\nBcCusWKeySDoX3+koJekWTf3yzhVdTDJVuBWYB1wdVXdl+SS4fkdwDuApwBXJQE4UFWbVq9sSdJy\nNFmzp6puBm4eO7ZjZPti4OJ2S5MktcU7aCWpgVlfxjHsJakDDHtJamBWZ/SHGPaStAyzGvqGvSR1\ngGEvSQ3M6oz+EMNekpZhVkPfsJekDjDsJamBWZ3RH2LYS9IyzGroG/aS1AGGvSQ1MKsz+kMMe0la\nhlkNfcNekjrAsJekBmZ1Rn+IYS9JyzCroW/YS1IHGPaS1MCszugPMewlaRlmNfQNe0nqAMNekhqY\n1Rn9IYa9JC3DrIb+xLBPspBkd5I9SS49ypj3Ds/fk+Ts9suUJK3EkmGfZB2wHVgAzgS2JDljbMz5\nwLOraiPwZuCqVap1qvX7/bUuYVXNc3/z3BvYX1tmdUZ/yKSZ/SZgb1U9UFUHgJ3A5rExrwY+DFBV\ndwAnJjm59UqnnP+gZtc89wb217ZZDf1JYX8KsG9kf//w2KQxp668NElSW9ZPON/0OSxNPu9Vr2p4\ntRn0la/A5z+/1lWsnnnub557A/try1e/OvhzVmf2qSUqT3IusFhVC8P9y4BHquqKkTEfAPpVtXO4\nvxs4r6oeGrvWjH6JJGltVdX4hHrZJs3s7wQ2JtkAPAhcAGwZG7ML2ArsHD45fHc86NsqVpL02CwZ\n9lV1MMlW4FZgHXB1Vd2X5JLh+R1VdVOS85PsBX4IXLjqVUuSlmXJZRxJ0nxY9Ttom9yUNe2SnJbk\nU0m+lOSLSf5oePypSW5L8tUkn0hy4sjnXDbseXeS31y76ptLsi7JXUluHO7PTX9JTkzy0ST3Jfly\nknPmpb9hrV9Kcm+SjyR5/Cz3luSaJA8luXfk2LL7SfKC4ddkT5K/PdZ9HM1R+vvr4ffmPUluSPLk\nkXPt9FdVq/bBYOlnL7ABOA64GzhjNR9zlfp4GvC84fYTgK8AZwB/Bbx1ePxS4N3D7TOHvR437H0v\n8Li17qNBn38K/AOwa7g/N/0xuBfkouH2euDJ89DfsL77gccP9/8R+P1Z7g14CXA2cO/IseX0c2jF\n4rPApuH2TcDCWve2RH+/cejvAXj3avS32jP7JjdlTb2q+npV3T3c/l/gPgb3Fxy+oWz4528PtzcD\n11fVgap6gMFf0KZjWvQyJTkVOB/4Ox59K+1c9DecJb2kqq6BwWtRVfU95qO/7wMHgBOSrAdOYPBm\nipntrapuBx4eO7ycfs5J8nTgiVX12eG4vx/5nDV1pP6q6raqemS4eweP3qvUWn+rHfZNbsqaKcN3\nJp3N4C/k5Hr0nUcPAYfuHH4Gg14PmYW+/wb4M+CRkWPz0t/pwDeTXJvkP5J8KMnPMwf9VdV3gCuB\n/2YQ8t+tqtuYg97GLLef8eNfYzb6BLiIwUwdWuxvtcN+rl79TfIE4GPAH1fVD0bP1eBnqaX6ndqv\nRZLfAr5RVXfxszfIAbPdH4Nlm+cD76+q5zN419jbRgfMan9JfgH4EwY/4j8DeEKS14+OmdXejqZB\nPzMryZ8DP62qj7R97dUO+68Bp43sn8b/fzaaGUmOYxD011XVx4eHH0rytOH5pwPfGB4f7/vU4bFp\n9avAq5P8J3A98GtJrmN++tsP7K+qzw33P8og/L8+B/39MvBvVfXtqjoI3AD8CvPR26jlfC/uHx4/\ndez4VPeZ5I0MllJfN3K4tf5WO+wP35SV5HgGN2XtWuXHbF2SAFcDX66q94yc2sXgxTCGf3585Phr\nkxyf5HRgI4MXU6ZSVb29qk6rqtOB1wKfrKo3MD/9fR3Yl+QXh4deDnwJuJHZ7283cG6Snxt+n74c\n+DLz0duoZX0vDv/Ovz9811WAN4x8ztRJssBgGXVzVf1k5FR7/R2DV55fyeDdK3uBy9byVfAV9PBi\nBmvZdwN3DT8WgKcC/wx8FfgEcOLI57x92PNu4BVr3cMyej2PR9+NMzf9Ac8FPgfcw2D2++R56Q94\nK4Mnr3sZvHh53Cz3xuCnyweBnzJ4ze/Cx9IP8ILh12Qv8N617muJ/i4C9gD/NZIv72+7P2+qkqQO\n8L8llKQOMOwlqQMMe0nqAMNekjrAsJekDjDsJakDDHtJ6gDDXpI64P8AFhkihkCRx6IAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0724c07080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(runoff[1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f0724902898>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEZRJREFUeJzt3X2MZXV9x/H3h12wKg0bYrM8zNo1BRJITF0wyxa13Bpp\nl01dNLEVU6XFpBoTorVN5UETpn9VbJpaakTSglltBRu0ZBGMrJZrMU1Wgd0VYVfYVNrdURYjggox\nLvLtH/ewXseZfTpn5s69+34lkzkPv3N/5zszuZ8533PvTKoKSdKx7bhRn4AkafQMA0mSYSBJMgwk\nSRgGkiQMA0kSLcMgya8l2Zpke5KHk/ztPOOuT/Jokh1J1rSZU5LUveVtDq6qnyb5vap6Nsly4GtJ\nXltVX3thTJINwBlVdWaS84EbgHXtTluS1KXWbaKqerZZPAFYBjw5a8hGYFMzdiuwIsnKtvNKkrrT\nOgySHJdkO7APuKeqHp415HRgz9D6XmCq7bySpO50cWXwfFW9isET/O8m6c0xLLMPazuvJKk7re4Z\nDKuqp5PcCbwa6A/tmgFWDa1PNdt+SRIDQpKOQlXN/oX7iLV9NdHLkqxoll8MXARsmzVsM3BZM2Yd\n8FRV7Zvr8apqYj+uvfbakZ+DtVmf9S3cx9RUAYtfX1faXhmcCmxKchyDYPl0VX0lybsBqurGqror\nyYYku4FngMtbzilJ6ljbl5Y+CJw7x/YbZ61f0WYeSdLC8h3Ii6TX6436FBbMJNcG1jfuJr2+rqTL\nnlMbSWqpnIskHalVq2DvXljsp7Ek1KhvIEuSJoNhIEkyDCRJhoEkCcNAkoRhIEnCMJAkYRhIkjAM\nJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJtAyDJKuS3JPkoSTfSvLeOcb0kjydZFvz\n8aE2c0qSure85fH7gfdX1fYkJwL3J9lSVTtnjftqVW1sOZckLVnj/l97W10ZVNXjVbW9Wf4JsBM4\nbY6hrf8/pyRp4XR2zyDJamANsHXWrgIuSLIjyV1JzulqTklaKjLmv/K2bRMB0LSIbgPe11whDHsA\nWFVVzya5GLgdOKuLeSVpqRj3NlHrMEhyPPA54F+r6vbZ+6vqx0PLX0zy8SQnV9WTs8dOT08fWO71\nevR6vbanJ0kTpd/v0+/3O3/cVIs4SxJgE/CDqnr/PGNWAk9UVSVZC/x7Va2eY1y1ORdJGqVVq2Dv\n3sW/QkhCVbVuUrW9MngN8Hbgm0m2NduuAV4OUFU3Am8B3pPkOeBZ4NKWc0rSkjPuv8u2ujLoklcG\nksbZ1BTMzIzvlYHvQJakDoz7q4kMA0nqwLg3NgwDSZJhIEldsE0kSRp7hoEkyTCQJBkGkiQMA0kS\nhoEkCcNAkoRhIEnCMJAkYRhIkjAMJKkT/qE6SdLYMwwkqQP+oTpJkm0iSdL4axUGSVYluSfJQ0m+\nleS984y7PsmjSXYkWdNmTklaisa9TbS85fH7gfdX1fYkJwL3J9lSVTtfGJBkA3BGVZ2Z5HzgBmBd\ny3klaUk5pttEVfV4VW1vln8C7AROmzVsI7CpGbMVWJFkZZt5JUnd6uyeQZLVwBpg66xdpwN7htb3\nAlNdzStJS8Gx3iYCoGkR3Qa8r7lC+JUhs9bnvKCanp4+sNzr9ej1el2cniQtuMVqE/X7ffr9fueP\nm2pZQZLjgS8AX6yqj86x/xNAv6pubdZ3ARdW1b5Z46rtuUjSqExNwczM4t87SEJVtb4uaftqogA3\nAQ/PFQSNzcBlzfh1wFOzg0CSxt2x3iZ6DfB24JtJtjXbrgFeDlBVN1bVXUk2JNkNPANc3nJOSVpy\nxr2x0bpN1BXbRJLG2THdJpIkDYx7m8gwkKQOjHtjwzCQJBkGktQF20SSJNtEkqTxZxhIUgdsE0mS\nbBNJksafYSBJMgwkqQveM5Akec9AkjT+DANJ6oBtIkmSbSJJ0vgzDCSpA7aJJEm2iSRJ4691GCS5\nOcm+JA/Os7+X5Okk25qPD7WdU5KWmnFvEy3v4DE+CfwT8KmDjPlqVW3sYC5JWpKO+TZRVd0L/PAQ\nw8Y8MyVpsi3GPYMCLkiyI8ldSc5ZhDklaVHZJjq0B4BVVfVskouB24Gz5ho4PT19YLnX69Hr9Rbh\n9CSpvcVqE/X7ffr9fuePm+qggiSrgTuq6pWHMfY7wHlV9eSs7dXFuUjSKExNwczM4t87SEJVtb4u\nWfA2UZKVyeACKslaBgH05CEOk6Sxcsy3iZLcAlwIvCzJHuBa4HiAqroReAvwniTPAc8Cl7adU5KW\nmnFvbHTSJuqCbSJJ48w2kSRp7NtEhoEkdWDcGxuGgSTJMJAkGQaSJAwDSRKGgSQJw0CShGEgScIw\nkCRhGEiSMAwkSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkuggDJLcnGRfkgcPMub6JI8m2ZFk\nTds5JWmp8d9ewieB9fPtTLIBOKOqzgTeBdzQwZyStKQc8//2sqruBX54kCEbgU3N2K3AiiQr284r\nSerOYtwzOB3YM7S+F5hahHkladGMe5to+SLNM/vLNOcF1fT09IHlXq9Hr9dbuDOSpA4tVpuo3+/T\n7/c7f9xUBxUkWQ3cUVWvnGPfJ4B+Vd3arO8CLqyqfbPGVRfnIkmjMDUFMzOLf+8gCVXV+rpkMdpE\nm4HLAJKsA56aHQSSpNFq3SZKcgtwIfCyJHuAa4HjAarqxqq6K8mGJLuBZ4DL284pSepWJ22iLtgm\nkjTObBNJksaeYSBJMgwkSYaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgSQJw0CShGEg\nScIwkCRhGEiSMAwkSRgGkiQMA0kSHYRBkvVJdiV5NMmVc+zvJXk6ybbm40Nt55QkdWt5m4OTLAM+\nBrwBmAG+kWRzVe2cNfSrVbWxzVySpIXT9spgLbC7qh6rqv3ArcAlc4xLy3kkSQuobRicDuwZWt/b\nbBtWwAVJdiS5K8k5LeeUJHWsVZuIwRP9oTwArKqqZ5NcDNwOnDXXwOnp6QPLvV6PXq/X8vQkabL0\n+336/X7nj5uqw3k+n+fgZB0wXVXrm/Wrgeer6rqDHPMd4LyqenLW9mpzLpI0SlNTMDMDi/00loSq\nat2Kb9smug84M8nqJCcAbwU2Dw9IsjJJmuW1DALoyV99KEnSqLRqE1XVc0muAL4ELANuqqqdSd7d\n7L8ReAvwniTPAc8Cl7Y8Z0lSx1q1ibpkm0jSODvW20SSpAlgGEiSDANJkmEgScIwkCRhGEiSMAwk\nSRgGkiQMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJJEB2GQZH2SXUke\nTXLlPGOub/bvSLKm7ZySpG61CoMky4CPAeuBc4C3JTl71pgNwBlVdSbwLuCGNnNKkrrX9spgLbC7\nqh6rqv3ArcAls8ZsBDYBVNVWYEWSlS3nlSR1qG0YnA7sGVrf22w71JiplvNKkjq0vOXxdZjjcjjH\nTU9PH1ju9Xr0er2jOilJmlT9fp9+v9/546bqcJ/P5zg4WQdMV9X6Zv1q4Pmqum5ozCeAflXd2qzv\nAi6sqn2zHqvanIskjdLUFMzMwGI/jSWhqmb/wn3E2raJ7gPOTLI6yQnAW4HNs8ZsBi6DA+Hx1Owg\nkCSNVqs2UVU9l+QK4EvAMuCmqtqZ5N3N/hur6q4kG5LsBp4BLm991pKkTrVqE3XJNpGkcXast4kk\nSRPAMJAkGQaSJMNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIw\nDCRJGAaSJAwDSRIt/gdykpOBzwK/CTwG/HFVPTXHuMeAHwE/B/ZX1dqjnVOStDDaXBlcBWypqrOA\nrzTrcymgV1VrDAJJWprahMFGYFOzvAl400HGtv5nzZKkhdMmDFZW1b5meR+wcp5xBXw5yX1J/rzF\nfJKkBXLQewZJtgCnzLHrg8MrVVVJap6HeU1VfS/JbwBbkuyqqnuP7nQlSQvhoGFQVRfNty/JviSn\nVNXjSU4FnpjnMb7XfP5+kv8A1gJzhsH09PSB5V6vR6/XO9T5S9Ixpd/v0+/3O3/cVM33C/0hDkw+\nAvygqq5LchWwoqqumjXmJcCyqvpxkpcCdwN/U1V3z/F4dbTnIkmjNjUFMzOw2E9jSaiq1vdl29wz\n+DBwUZJHgNc36yQ5LcmdzZhTgHuTbAe2Al+YKwgkSaN11FcGXfPKQNI4O5avDCRJE8IwkKQOjHtj\nwzCQJBkGktSFjPnfWTAMJKkDtokkSWPPMJCkDtgmkiTZJpIkjT/DQJI6YJtIkmSbSJI0/gwDSeqA\nbSJJ0tgzDCRJhoEkyTCQJGEYSJIwDCRJtAiDJH+U5KEkP09y7kHGrU+yK8mjSa482vkkSQunzZXB\ng8Cbgf+ab0CSZcDHgPXAOcDbkpzdYs6x1e/3R30KC2aSawPrG3eTXl9XjjoMqmpXVT1yiGFrgd1V\n9VhV7QduBS452jnH2ST/QE5ybWB9427S6+vKQt8zOB3YM7S+t9kmSVpClh9sZ5ItwClz7Lqmqu44\njMc/oj/d9MY3Hsno8fLtb8P994/6LBbGJNcG1jfuFqu+mZmFn2MhpVr+qb0k9wB/VVUPzLFvHTBd\nVeub9auB56vqujnGjvnf/JOk0aiq1n8Z6aBXBkdgvhO5DzgzyWrgu8BbgbfNNbCLYiRJR6fNS0vf\nnGQPsA64M8kXm+2nJbkToKqeA64AvgQ8DHy2qna2P21JUpdat4kkSeNv5O9AnoQ3pSVZleSe5k14\n30ry3mb7yUm2JHkkyd1JVgwdc3VT864kvz+6sz88SZYl2ZbkjmZ9kmpbkeS2JDuTPJzk/Amr7+rm\nZ/PBJJ9J8qJxri/JzUn2JXlwaNsR15PkvOZr8miSf1zsOuYzT31/1/x87kjy+SQnDe3rpr6qGtkH\nsAzYDawGjge2A2eP8pyOso5TgFc1yycC3wbOBj4CfKDZfiXw4Wb5nKbW45vadwPHjbqOQ9T4l8C/\nAZub9UmqbRPwzmZ5OXDSpNTXnOP/AC9q1j8L/Ok41we8DlgDPDi07UjqeaEj8nVgbbN8F7B+1LUd\npL6LXvg+AB9eiPpGfWUwEW9Kq6rHq2p7s/wTYCeD91NsZPBEQ/P5Tc3yJcAtVbW/qh5j8A1cu6gn\nfQSSTAEbgH/hFy8WmJTaTgJeV1U3w+A+V1U9zYTUB/wI2A+8JMly4CUMXswxtvVV1b3AD2dtPpJ6\nzk9yKvDrVfX1Ztynho4Zqbnqq6otVfV8s7oVmGqWO6tv1GEwcW9Ka145tYbBN2xlVe1rdu0DVjbL\npzGo9QVLve5/AP4aeH5o26TU9grg+0k+meSBJP+c5KVMSH1V9STw98D/MQiBp6pqCxNS35AjrWf2\n9hnGo06AdzL4TR86rG/UYTBRd6+TnAh8DnhfVf14eF8NrtUOVu+S/Fok+UPgiaraxjwvIR7X2hrL\ngXOBj1fVucAzwFXDA8a5viS/BfwFgxbCacCJSd4+PGac65vLYdQztpJ8EPhZVX2m68cedRjMAKuG\n1lfxy2k2NpIczyAIPl1Vtzeb9yU5pdl/KvBEs3123VPNtqXoAmBjku8AtwCvT/JpJqM2GPy87a2q\nbzTrtzEIh8cnpL5XA/9dVT+owUu9Pw/8DpNT3wuO5Odxb7N9atb2JV1nkj9j0K79k6HNndU36jA4\n8Ka0JCcweFPa5hGf0xFLEuAm4OGq+ujQrs0MbtbRfL59aPulSU5I8grgTAY3e5acqrqmqlZV1SuA\nS4H/rKp3MAG1weB+D7AnyVnNpjcADwF3MAH1AbuAdUle3PycvoHBe34mpb4XHNHPY/N9/1HzyrEA\n7xg6ZslJsp5Bq/aSqvrp0K7u6lsCd84vZvDqm93A1aM+n6Os4bUM+unbgW3Nx3rgZODLwCPA3cCK\noWOuaWreBfzBqGs4zDov5BevJpqY2oDfBr4B7GDwm/NJE1bfBxgE3IMMbq4eP871MbhC/S7wMwb3\nHC8/mnqA85qvyW7g+lHXdZD63gk8Cvzv0PPLx7uuzzedSZJG3iaSJC0BhoEkyTCQJBkGkiQMA0kS\nhoEkCcNAkoRhIEkC/h+1qd5XR1P7lwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f07249335f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(run_depth[1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f07226f70f0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEACAYAAABPiSrXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEydJREFUeJzt3H+wHeV93/H3xwg8pGbATDuIHzJWG3mCapoaYuTYTX0a\n21TpdCwy0xSYhjo2bcfVtE5/jG1kpkV/ucRpk+C28EcDwbiFhHFcBk+wajnlpHVaIycBW1jISNMK\nkGyw49aQtE0slW//OHvk49srsffq7t67575fM2e0++ye3ecrLvqc59nz3FQVkiS18arV7oAkaTgM\nDUlSa4aGJKk1Q0OS1JqhIUlqzdCQJLXWWWgk2ZTk0SRfTfJkkg807Rcm2Zvk6SSfS3LBzHt2JTmU\n5GCSa7vqmyRpedLVOo0kG4GNVfVEktcAvwtcB7wX+P2q+liSDwOvrapbkmwF7gfeDFwKfB54Q1W9\n3EkHJUlL1tlIo6qer6onmu0/BJ5iEgbvBj7RnPYJJkECsAN4oKqOV9UR4DBwTVf9kyQtXS/PNJK8\nHngT8BhwUVW90Bx6Abio2b4EODrztqNMQkaStEZ0HhrN1NSvAz9bVX8we6wmc2Onmx/zd5xI0hqy\nocuLJzmbSWB8sqoeappfSLKxqp5PcjHwzab9GLBp5u2XNW2z1zNEJGkZqiorcZ0uvz0V4G7gQFX9\n0syhh4H3NNvvAR6aab8hyTlJNgNbgH0Lr1tVc/u67bbbVr0P1md966229VDfSupypPE24KeBryR5\nvGnbBdwOPJjkZuAI8NcBqupAkgeBA8AJYGetdLWSpDPSWWhU1Rc49Ujmnad4z0eBj3bVJ0nSmXFF\n+BoyGo1Wuwudsr7hmufaYP7rW0mdLe7rQhJnrCRpiZJQa/1BuCRp/hgakqTWDA1JUmuGhiSpNUND\nktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIkloz\nNCRJrRkakqTWDA1JUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTW1lRoJNme5GCS\nQ0k+vNr9kSR9vzUTGknOAv4VsB3YCtyY5IrV7ZUkrayvfW21e3Bm1kxoANcAh6vqSFUdB34V2LHK\nfZKkFfVDPwTf+MZq92L51lJoXAo8N7N/tGmTpLly4sRq92D51lJo1Gp3QJJ0ehtWuwMzjgGbZvY3\nMRltfJ/du3ef3B6NRoxGo677JUmDMh6PGY/HnVw7VWvjA36SDcDXgHcAXwf2ATdW1VMz59Ra6a8k\nLUcCzz4Lmza98rkrd89QVVmJa62ZkUZVnUjy94D/AJwF3D0bGJKk1bdmRhptONKQNHRDH2mspQfh\nkqQ1ztCQJLVmaEiSWjM0JEmtGRqSpNYMDUlSa4aGJKk1Q0OS1JqhIUlqzdCQJLVmaEiSWjM0JEmt\nGRqSpNYMDUlSa4aGJKk1Q0OS1JqhIUlqzdCQJLVmaEiSWjM0JEmtGRqSpNYMDUlSa4aGJKk1Q0OS\n1JqhIUlqzdCQJLVmaEiSWjM0JEmtGRqSpNY6CY0kP5/kqSRfTvLpJOfPHNuV5FCSg0munWm/Osn+\n5tgdXfRLknRmuhppfA74s1X1w8DTwC6AJFuB64GtwHbgziRp3nMXcHNVbQG2JNneUd8kScvUSWhU\n1d6qernZfQy4rNneATxQVcer6ghwGNiW5GLgvKra15x3H3BdF32TJC1fH8803gc80mxfAhydOXYU\nuHSR9mNNuyTNjarV7sGZ27DcNybZC2xc5NBHquozzTm3At+tqvuXe5+Fdu/efXJ7NBoxGo1W6tKS\n1Iuuw2M8HjMejzu5dqqj3if5GeBvA++oqj9q2m4BqKrbm/09wG3AM8CjVXVF034j8Paqev+Ca1ZX\n/ZWkrlXBq14FzzwDr3tdf/dNQlXllc98ZV19e2o78EFgxzQwGg8DNyQ5J8lmYAuwr6qeB15Ksq15\nMH4T8FAXfZOk1TL9zDvkz77Lnp56Bf8SOAfY23w56r9W1c6qOpDkQeAAcALYOTN02AncC5wLPFJV\nezrqmyRpmTqbnuqC01OShuzll+Gss+DIEbj88v7uu+anpyRJ/795mJ4yNCRJrRkaktSTIY8wpgwN\nSerZkMPD0JAktWZoSFJPfBAuSVpXDA1JUmuGhiT1xOkpSdK6YmhIUk+GPMKYMjQkqWdDDg9DQ5LU\nmqEhST3xQbgkaV0xNCRJrRkaktQTp6ckSeuKoSFJPRnyCGPK0JCkng05PAwNSVJrhoYk9cQH4ZKk\ndcXQkCS1ZmhIUk+cnpIkrSuGhiT1xJGGJGld6TQ0kvzjJC8nuXCmbVeSQ0kOJrl2pv3qJPubY3d0\n2S9J0vJ0FhpJNgHvAp6ZadsKXA9sBbYDdyZJc/gu4Oaq2gJsSbK9q75J0mpweur0fgH40IK2HcAD\nVXW8qo4Ah4FtSS4Gzquqfc159wHXddg3SdIydBIaSXYAR6vqKwsOXQIcndk/Cly6SPuxpl2S5saQ\nRxhTG5b7xiR7gY2LHLoV2AVcO3v6cu+z0O7du09uj0YjRqPRSl1aknrRdXiMx2PG43En106tcO+T\nvBH4TeB/N02XMRk5bAPeC1BVtzfn7gFuY/Lc49GquqJpvxF4e1W9f8G1a6X7K0l9efFFuOAC2L8f\n3vjG/u6bhKpakQ/vKz49VVVPVtVFVbW5qjYzmXa6qqpeAB4GbkhyTpLNwBZgX1U9D7yUZFvzYPwm\n4KGV7pskraZ5eBC+7OmpJTj511NVB5I8CBwATgA7Z4YOO4F7gXOBR6pqTw99kyQtwYpPT3XJ6SlJ\nQ/ad78BrXwtf+QpceWV/913T01OSpMXNw/SUoSFJas3QkKSeDHmEMWVoSFLPhhwehoYkqTVDQ5J6\n4oNwSdK6YmhIklozNCSpJ05PSZLWFUNDknoy5BHGlKEhST0bcngYGpKk1gwNSeqJD8IlSeuKoSFJ\nas3QkKSeOD0lSVpXDA1J6smQRxhThoYk9WzI4WFoSJJaMzQkqSc+CJckrSuGhiSpNUNDknri9JQk\naV0xNCSpJ0MeYUwZGpLUsyGHR2ehkeTvJ3kqyZNJfm6mfVeSQ0kOJrl2pv3qJPubY3d01S9J0vJt\n6OKiSf4S8G7gz1XV8SR/qmnfClwPbAUuBT6fZEtVFXAXcHNV7UvySJLtVbWni/5J0mrwQfip/V3g\nn1XVcYCq+lbTvgN4oKqOV9UR4DCwLcnFwHlVta857z7guo76Jklapq5CYwvwF5N8Mck4yY807ZcA\nR2fOO8pkxLGw/VjTLklaQ5Y9PZVkL7BxkUO3Ntd9bVW9JcmbgQeBP73ce83avXv3ye3RaMRoNFqJ\ny0pS5/qanhqPx4zH406uneqg90k+C9xeVb/V7B8G3gL8LYCqur1p3wPcBjwDPFpVVzTtNwJvr6r3\nL7huddFfSerDc8/B614Hv/3b8Na39nffJFRVVuJaXU1PPQT8OECSNwDnVNXvAw8DNyQ5J8lmJtNY\n+6rqeeClJNuSBLipuYYkzY15+MzbybengHuAe5LsB74L/E2AqjqQ5EHgAHAC2DkzdNgJ3AucCzzi\nN6ckzashh0cn01NdcXpK0pA9+yxcfjl84Qvwtrf1d98hTE9JkhZwnYYkaV0xNCSpJ0MeYUwZGpLU\nsyGHh6EhSWrN0JCkngx5hDFlaEhSz4YcHoaGJKk1Q0OSeuI6DUnSumJoSFJPhjzCmDI0JKlnQw4P\nQ0OS1JqhIUk9GfIIY8rQkKSeDTk8DA1JUmuGhiT1xHUakqR1xdCQpJ4MeYQxZWhIUs+GHB6GhiSp\nNUNDknoy5BHGlKEhST0bcngYGpKk1gwNSeqJ6zQkSeuKoSFJPRnyCGPK0JCkng05PDoJjSTXJNmX\n5PEkX0ry5plju5IcSnIwybUz7Vcn2d8cu6OLfkmSzkxXI42PAf+kqt4E/NNmnyRbgeuBrcB24M4k\nad5zF3BzVW0BtiTZ3lHfJGlVDHmEMdVVaHwDOL/ZvgA41mzvAB6oquNVdQQ4DGxLcjFwXlXta867\nD7iuo75J0qoacnhs6Oi6twBfSPLPmQTTjzbtlwBfnDnvKHApcLzZnjrWtEuS1pBlh0aSvcDGRQ7d\nCnwA+EBV/fskPwXcA7xrufeSpHkwD+s0lh0aVXXKEEjyb6vqnc3up4BfbraPAZtmTr2MyQjjWLM9\n236MRezevfvk9mg0YjQaLbHnkjTfxuMx4/G4k2unOoi8JL8H/MOq+q0k7wBur6o3Nw/C7weuYTL9\n9HngB6uqkjzGZISyD/gN4ONVtWfBdauL/kpSH558Eq68Ej77Wdje41d9klBVeeUzX1lXzzT+DvCv\nk7wa+D/NPlV1IMmDwAHgBLBzJgV2AvcC5wKPLAwMSZoXQ/7s20loVNXvANtOceyjwEcXaf9d4Mou\n+iNJWhmuCJekngx5hDFlaEhSz4YcHoaGJPVkyGExZWhIUs+GHB6GhiSpNUNDknoy5BHGlKEhST0b\ncngYGpKk1gwNSerJkEcYU4aGJPVsyOFhaEhST4YcFlOGhiT1bMjhYWhIklozNCSpJ0MeYUwZGpLU\nsyGHh6EhSWrN0JCkngx5hDFlaEhSz4YcHoaGJPVkyGExZWhIUs+GHB6GhiSpNUNDknoy5BHGlKEh\nST0bcngYGpKk1gwNSerJkEcYU4aGJPVsyOFhaEhST4YcFlOGhiT1bMjhsezQSPJTSb6a5P8muWrB\nsV1JDiU5mOTamfark+xvjt0x0/7qJL/WtH8xyeXL7ZckqTtnMtLYD/wk8J9mG5NsBa4HtgLbgTuT\npDl8F3BzVW0BtiTZ3rTfDHy7af9F4OfOoF+DNR6PV7sLnbK+4Zrn2qC/+oY8wphadmhU1cGqenqR\nQzuAB6rqeFUdAQ4D25JcDJxXVfua8+4Drmu23w18otn+deAdy+3XkPk/5rDNc33zXBv0X9+Qw6OL\nZxqXAEdn9o8Cly7Sfqxpp/nzOYCqOgG8mOTCDvomSToDG053MMleYOMihz5SVZ/ppkunt2nTaty1\nHy++CHffvdq96I71Ddc81wb91ffHfzz5c8gjDarqjF7Ao8BVM/u3ALfM7O8BtjEJn6dm2m8E7po5\n5y3N9gbgW6e4V/ny5cuXr6W/zvTf+unrtCONJcjM9sPA/Ul+gcm00xZgX1VVkpeSbAP2ATcBH595\nz3uALwJ/DfjNxW5SVVmsXZLUj2WHRpKfZPKP/p8EfiPJ41X1E1V1IMmDwAHgBLCz6uRgbCdwL3Au\n8EhV7Wna7wY+meQQ8G3ghuX2S5LUnXzv33NJkk5vMCvCk2xvFgseSvLh1e7PUiXZlOTRZkHkk0k+\n0LRfmGRvkqeTfC7JBTPvWXSR5FqW5Kwkjyf5TLM/N/UluSDJp5I8leRAkm1zVt+u5udzf5L7m0W3\ng6wvyT1JXkiyf6ZtybWcakHyajtFfT/f/Gx+Ocmnk5w/c2zl6luphyNdvoCzmKz3eD1wNvAEcMVq\n92uJNWwE/nyz/Rrga8AVwMeADzXtHwZub7a3NnWe3dR9GHjVatfRos5/BPw74OFmf27qY7KW6H3N\n9gbg/Hmpr+njfwNe3ez/GpPnjIOsD/gx4E3A/pm2pdQynYXZB1zTbD8CbF/t2k5T37um/w2A27uq\nbygjjWuAw1V1pKqOA7/KZBHhYFTV81X1RLP9h8BTTL4oMLuw8RN8b8HjYoskr+m100uU5DLgrwC/\nzPe+HDEX9TWf2n6squ4BqKoTVfUic1If8BJwHPiBJBuAHwC+zkDrq6r/DPzPBc1LqeWVFiSvqsXq\nq6q9VfVys/sYcFmzvaL1DSU0Ti7+a0wXDA5Sktcz+ZTwGHBRVb3QHHoBuKjZPtUiybXsF4EPAi/P\ntM1LfZuBbyX5lSS/l+TfJPkTzEl9VfU/gH8BPMskLL5TVXuZk/oaS63ldAuS17r3MRk5wArXN5TQ\nmJun9Ulew+RXpfxsVf3B7LGajBFPV+ua/XtI8leBb1bV43z/V7BPGnJ9TKajrgLurKqrgP/FZE3S\nSUOuL8mfAf4Bk+mLS4DXJPnp2XOGXN9CLWoZrCS3At+tqvu7uP5QQuMYMLsWfBPfn5CDkORsJoHx\nyap6qGl+IcnG5vjFwDeb9oU1X9a0rVVvBd6d5L8DDwA/nuSTzE99R4GjVfWlZv9TTELk+Tmp70eA\n/1JV367Jr/L5NPCjzE99sLSfxaNN+2UL2td0jUl+hskU8d+YaV7R+oYSGr/D5Lfivj7JOUx+i+7D\nq9ynJUkSJutRDlTVL80cmi5spPnzoZn2G5Kck2QzzSLJvvq7VFX1karaVFWbmayz+Y9VdRPzU9/z\nwHNJ3tA0vRP4KvAZ5qA+4CDwliTnNj+r72Sy1mpe6oMl/iw2/81far4lFyYLkh9aeNG1IpPfGv5B\nYEdV/dHMoZWtb7W/BbCEbwv8BJNvHB0Gdq12f5bR/7/AZK7/CeDx5rUduBD4PPA08Dnggpn3fKSp\n9yDwl1e7hiXU+na+9+2puakP+GHgS8CXmXwSP3/O6vsQkyDcz+RB8dlDrY/JaPfrwHeZPA9973Jq\nAa5u/j4OAx9f7bpOU9/7gEPAMzP/vtzZRX0u7pMktTaU6SlJ0hpgaEiSWjM0JEmtGRqSpNYMDUlS\na4aGJKk1Q0OS1JqhIUlq7f8BrL1bFJajQ5AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0724914d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(run_temp[1,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create netcdf file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nemo = nc.Dataset('../rivers/river_cnst.nc', 'w')\n",
    "nemo.description = 'Constant Yearly Average, One River' \n",
    "# dimensions\n",
    "nemo.createDimension('x', xmax)\n",
    "nemo.createDimension('y', ymax)\n",
    "nemo.createDimension('time_counter', None)\n",
    "    \n",
    "# variables\n",
    "y = nemo.createVariable('y','float32',('y','x'),zlib=True)\n",
    "y = Y\n",
    "x = nemo.createVariable('x','float32',('y','x'),zlib=True)\n",
    "x = X\n",
    "    # time\n",
    "time_counter = nemo.createVariable('time_counter', 'float32', ('time_counter'),zlib=True)\n",
    "time_counter.units = 'non-dim'\n",
    "time_counter[0] = 1\n",
    "# runoff\n",
    "rorunoff = nemo.createVariable('rorunoff', 'float32', ('time_counter','y','x'), zlib=True)\n",
    "rorunoff._Fillvalue = 0.\n",
    "rorunoff._missing_value = 0.\n",
    "rorunoff._units = 'kg m-2 s-1'\n",
    "rorunoff[0,:] = runoff\n",
    "# depth\n",
    "rodepth = nemo.createVariable('rodepth','float32',('y','x'),zlib=True)\n",
    "rodepth._Fillvalue = -1.\n",
    "rodepth.missing_value = -1.\n",
    "rodepth.units = 'm'\n",
    "rodepth[:] = run_depth\n",
    "# temp\n",
    "rotemper = nemo.createVariable('rotemper','float32',('y','x'),zlib=True)\n",
    "rotemper._Fillvalue = -999.\n",
    "rotemper.missing_value = -999.\n",
    "rotemper.units = 'deg C'\n",
    "rotemper[:] = run_temp\n",
    "nemo.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}