smb_downscaling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CLM-CISM downscaling, 20 Oct 2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CLM forcing fields (acab or SMB, and temperature) are supplied at a coarse horizontal resolution and at 10 vertical levels or elevation classes.\n",
    "\n",
    "Here I explore some possibilities regarding horizontal and vertical interpolation of these fields to the CISM grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two cells are selected along the Greenland K-transect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "coords={} # Coordinates (ilat,ilon) of CLM gridcells\n",
    "coords['cellA']=(167,248)\n",
    "coords['cellB']=(167,249)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some test data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yearly average ablation (AB) and accumulation (AC) are taken from an FG run (casename = cam6clm5_beta07_col_119rest_lwlap) with output at each elevation class.\n",
    "\n",
    "One year averaged data for now, just for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Elevation class midpoints\n",
    "elev_mec = np.asarray([100, 300, 550, 850, 1150, 1450, 1800, 2250, 2750, 3500])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ab={}\n",
    "ab['cellA']=np.asarray([  1.04427941e+01 ,  7.75391426e+00 , 5.62771183e+00,   2.39465291e+00,\n",
    "   7.41343991e-01 ,  7.64908130e-02 ,  5.09881154e-02  , 1.36675815e-02,\n",
    "  -1.02899391e-02 , -3.26783274e-02]) * 365 / 1000. # convert from mm/day to m/year\n",
    "ab['cellB']=np.asarray([  1.53967969e+01,   1.20004906e+01,   8.89361040e+00 ,  3.97812184e+00,\n",
    "   1.51039940e+00  , 1.65159715e-01 ,  1.29554068e-01  , 7.93311018e-02,\n",
    "   3.70435832e-02,  -6.42650455e-03]) * 365 / 1000. # convert from mm/day to m/year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ac={}\n",
    "ac['cellA']=np.asarray([ 0.48337985,  0.62594054,  0.77821437,  0.9810329 ,  1.04418045,  1.12301125,\n",
    "  1.14667178 , 1.15020953 , 1.15020964,  1.15020964]) * 365 / 1000. # convert from mm/day to m/year\n",
    "ac['cellB']=np.asarray([ 0.5283285 ,  0.65806491 , 0.80689576 , 1.01173176 , 1.11583912,  1.19848984,\n",
    "  1.24245121,  1.25205707 , 1.25205707 , 1.25205707]) * 365 / 1000. # convert from mm/day to m/year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SMB can be calculated as AB-AC:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "smb = {}\n",
    "smb['cellA'] = ac['cellA']-ab['cellA']\n",
    "smb['cellB'] = ac['cellB']-ab['cellB']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10c757310>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEZCAYAAABmTgnDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXJ2wRFQWq4gYuGJfiWkXFIqmtS/0qaimt\n+1IV64LYKj9xaYnfakXcra2KuIDfKrUuKNUiuATBBVFAUdHgBoiAAgKyBEjy+f1xbmAg2ySZmTuT\nvJ+Pxzwyc+becz5zSebDPffcc8zdERERqU1e3AGIiEj2U7IQEZE6KVmIiEidlCxERKROShYiIlIn\nJQsREamTkoVkDTPb1sxeN7NlZnZr3PHEwczOMbOJDdy3l5nNTXVMIqBkIY1gZj81szfMbKmZLTKz\niWb2k0ZU2Q/41t23cveBqYozBzXm5ifdOCVp0TLuACQ3mdmWwBjgIuDfQGugJ7CmgfUZ0AX4OFUx\nikjq6MxCGqoAcHd/0oM17v6yu38IYGaDzeyxyo3NrIuZVZhZXvT6NTO70cwmmdlKYARwDnC1mS03\ns6PM7BAze9PMvjezeWb2NzNrmVDnj81snJktNrP5ZjYoKjczG2Rmn5nZd2Y2ysy2ru5DVHbdmNlA\nM1sYtXOSmf3SzD6NzpiuSdi+rpjujOpZZmbvm9k+UfnxZvZR9Nnmmtkfazm2eVG9S83sYzM7KqH+\nc6Oy5dHn61dTJWZ2dbTNcjP70MxOTnjvnOhM8FYzW2Jmn5vZcQnvtzezh6PPuNjMnkl47wQzmxYd\ng0lmtm8tn0WaCnfXQ496P4Atge+AR4HjgK03eX8wMDLhdRegHMiLXr8GfAXsRfhPSwvgEeB/E/Y5\nCOgOGNAZ+Ai4PHpvC+Ab4ArCWc3mwCHRewOAN4HtgVbAfcDjNXyOXsA64LoohguAb4H/A9oC+wCr\ngC5JxHQMMAXYMnq9J7Bd9PwboEf0fCvggBriOSeK5/Iont8ASyuPL/BLYJfoeU9gZWVd0WeZk1BX\nn4T2+wIrEl6fQzgL/F30WX4PzEvY9wXgCaBdFEfPqPxAYCFwcLTfWcCXQKu4fyf1SO8j9gD0yN1H\n9GX4MDAHWAs8B2wTvZdMsijapL6NkkU17Q0Ano6enwq8V8N2HwM/S3i9fRRfXjXb9oq+cC16vQVQ\nARycsM27QO8kYvoZ8AlwaGV9Cdt9BVxYmUhq+YznAF9vUjYZOKOG7Z8F+id8ljm11D0NODGhnZKE\n9zaLPve2QCegDGhXTR3/AG7YpOyTymSiR9N9qBtKGszdP3X337l7Z6AbsANwVz2qqHXkjpntYWZj\noi6mpcBNwI+it3cGPq9h1y7As1H3yhJC8lgHbFfD9os9+tYDVkc/v014fzUhidQak7u/BtwL/B1Y\naGb3m9kWUR19gP8BZkddcIfV8tHnbfJ6NuHYEnWPvRV1DX1PONP40aYVRNuendBd9D3w4022XVD5\nxN0rP/cWhGO7xN2XV1NtF+DKymMb1btTZXzSdClZSEq4ewmhS6pbVLSS0I1Tafvqdquj2vuAmcDu\n7r41oavIovfmArvXsN8c4Jfu3iF6tHf3zd19ft2fpE61xYS73+vuBxO6r/YEBkbl77n7ycA2hDOw\nJ2tpY8dNXncGvjGz1sBTwFDCGVx74L+J7Vcys87AMOCS6PO3J3SZVdm2GnOBDmbWrob3btrk2G7h\n7v9Kol7JYUoW0iBmtqeZ/dHMdoxe7wycBrwVbTIdONLMdjazrYBBDWhmS2C5u68ys72AixPe+w/Q\nycwuN7PWZraFmXWP3nsA+Gv0hYmZbWNmvRvQfr1iMrODzax7dMF7NVAKVJhZKzM73czauXs58AOh\nS64m25lZfzNraWZ9Cdd1XiBcm2kNLHL3CjP7JeE6SXU2J3QrLTKzPDM7jw2JvFbuvoCQhP5hZltH\ncfSM3n4Q+H3lsTazzaOL95snU7fkLiULaagfCH3zk83sB8IF5Q+AqwDc/WXgX1HZFMIw20TVnVVs\nWnYVcIaZLSckgFHrN3RfARwN9CZ0p5QAhdHbdxP+9z7OzJZFsXUneZvGkfi6xpgIF4MfBJYQLvou\nAipvLjwL+DLquuoHnF5L+28De0T7/wXo4+5Lo898OfDvqHvt1OhzVv0A7jOB26O6FhC6oCbV0uam\nn/MswnWLTwgXtAdE9b5HuPZybxRDCeH6hzRxlRf14gvA7CHgBGChu+8XazAiIlKtbDizeAQ4Nu4g\nRESkZrEnC3efBHwfdxwiIlKz2JOFiIhkPyULERGpU05MJGhmmklTRKQB3D2Ze2vqlC1nFkYdNwvF\nfat7Yx6DBw+OPYbmGn8ux67443/kevypFHuyMLPHCePgC8xsTnTzkIiIZJHYu6Hcvbabk0REJAvE\nfmbRHBQWFsYdQqPkcvy5HDso/rjlevypFPsd3MkwM8+FOEVEsomZ4U3sAreIiGQxJQsREamTkoWI\niNRJyUJEROqkZCEiInVSshARkTopWYiISJ2ULEREpE6xT/chIlIfQ/r1o7SkpEp5fkEBg4YNiyGi\n7FLT8Wks3cEt0swcv1d3Wi/4tkr52k7b8uIn78QQUf2ctv1uPLHgy6rlnXbliflfxBBRdkk8Pkbq\npijXmYXkvGz48suGGJLVesG3jF42u0r5yTHEgjtUVIRHeXl4VD6vrqy8nPIVa6qvakUpTJsWtk2s\nt/IRV1mG29520cK0/FMpWUjOy4Yvv5piOMUdVq2CsrKNH+vWVS2rpdzXlVGxZh0Va8vCY10ZvjaU\n+9pQ7usSHmVlsHbd+teVdVBWxk4rFlX7GXb5YSGLuh8fvpQ9fDlbRfiitvJyqCjHKiqin1F5RTnm\nCe97VJb4fmK5b1yW5xW4GRXWArc8PC/8rLAWVZ+TR0VeC3ZbUf2X4c4rFjHvuN/hWKiLvGj/vI3K\nKiwvem4bXkc/wTZ67YTtKjbab5My8qiooSyx3K1qWdi+ZSjzhLLE7WyTsmq2q0goW1g+FliVwt/u\nQMlCmqwWa9cx6b4ZeOkafHVp+Bk9WLMG1q6B0jXY2vBg7Rryoue2bg0t1paSt24NeWVraLFuDS3K\nokf5GlpW/ixfQ6vyNey9em61Mey7fA6rt/gR5daSMlqu/xkerdY/X0dCubdkLa0o85as8w3vldOS\nMmtFubWkPK8lFdZy/fNya0VFXksq8qL38lpSYfnhZ4vwnrcI5UsqngFWVon1G9+SO9ddtv5Lm5bR\nl3deC8gLX9q02PClXvl8o595LbAWG++3aXnlT1qE9/NaGGZh02R+Th+yC6ypmphntNmB//vjNMxI\n6gHJbZeq/cxCt1DLNLc55djrYeXixv8BbULJQnLW2rUwYQKUrqr+elbX1QvZaeBplLVoQ1mLNqxr\nmU95yzaUt2gTfrZsQ3mrNlS0bENF6zZUtGpDWas2+GbtqdiqDbRug7fJx1u3gTbhYfltIL8N1qYN\neZuF15bfhk8G9ILV86rEMHXzLnw66StaRN+NiY+8PGjTAtpWU17dtqly8tY3wbIlVY9nu7bcNO34\n1DWUJiffBVTTE5WfD1dfnfFwsk7LNH2rK1lITlm8GF58EcaMgXEvORds/wL7r6v6JQ3w6VY7scvS\nDzMS1z0DW8LqquUtW8IBB2QkhGbju806Urhslyrl6zb7IfPBZKGNj8+ElNWrZCFZzR0+/TQkh+ef\nhw8+gKOOgjN/MpPHFv2BNt98xVltO8Cq7+IONWes7bRttddz1nbaNuOxNMSPT+xLSUlp1fKC/Bii\nyT4bHZ8JqUsWGjorWWfdOpg0KSSIMWNg9Wo48UTo3Rt+duBS8m+5Af7v/+Daa+Gyyzh+3yNiH4mU\nS6OhpPlI5eJHShaSFb7/HsaODclh7FjYffeQIE48MXTjWEU5DB8OgwfDSSfBX/4C2+bG/4RF4pLK\nZKFuKInNZ59tOHt4913o1Sskh9tugx12SNhwwgQYMADatYP//hcOPDC2mEWaK51ZSEr16zek2v7k\ngoJ87rtvEG+9teH6w9KlcMIJoXvp5z+Htm032Wn2bBg4ECZPhltvhb59N4wjFJE66cxCslZJSSkT\nJhRVKZ85s4hnn4UddwzJYeRI+MlPahgSunIlDB0K994bzigefbSaTCIimRR7sjCz44C7CDPgPuTu\nt8QckqRBu3bwyivQuXMtG7nDqFFhsPwRR4SpG2rdQUQyJdZkYWZ5wL3Az4FvgClm9py7fxJnXNJw\nZWXVl++4Yx3f+++9F84iVq2Cf/4TevZMS3wi0jBxr2fRHZjl7rPdfR0wCjgp5pikgYqLYcqUeu60\ncCFccAH8z//AueeGCpQoRLJO3MliRyBxUp2vozLJIaWlcNVVcMYZUFCQ5E5r18Ltt0O3brD11uHO\nuwsuCHNbiEjWif2aRbKKiorWPy8sLKSwsDC2WGSD99+HM8+EPfcMz6+9Np+OHYuqbFeQeHftCy/A\nH/4Ae+wR7r7bc8/MBSzShBUXF1NcXJyWumMdOmtmhwFF7n5c9HoQ4Jte5NbQ2exTXh5ODG69Fe64\nIySMOke1fvIJ/PGP8PnncOedcHz2T1onksua0tDZKUBXM+sCzAdOBU6LNySpy1dfwdlnh+Tw7rvQ\npcuG96pd0rGsjPxlyxi0YAFccw2MHg2tW2c0ZhFpnFiThbuXm9llwDg2DJ2dGWdMUjP3cH/EVVeF\n0a1/+EPVSwylJSUUVTN5WVGnTvDRR5qiQyRHxX1mgbuPBdRpneUWLYKLLoJZs8L9EvvtV88K9txT\niUIkh8U9GkpywIsvwv77w267wTvvNCBRiEjOi/3MQrLXypWhy+m//w33ySU1AG3FinSHJSIx0JmF\nVGvy5DC568qVYUhsUonik0/C6kQi0uTozEI2sm4d3HQT3HdfmMevb98kd/ziCzj6aPJ79KComrfz\nk75bT0SykaYol/U+/RTOOgs6dICHH95kTYnazJ0LRx4ZphO/5JK0xigiyUvlfRbqhhLc4R//CBO9\nnnNOuEaRdKJYsCAsRnHppUoUIk2YuqGaufnz4Xe/C0NjJ02Cvfaqx86LFsHRR4fbt6+6Km0xikj8\ndGbRjD39dLiI3b07vPlmPRPF0qVw7LFhyo4//SltMYpIdtCZRTO0bBlcfnlIEKNHw2GH1bOCFStC\nkujRA4YM0VKnIs2AkkUTVt162EuXwuef53PGGYOYPh0237yela5eHdZF3XtvuPtuJQqRZkLJogmr\naT3sbt2KuP/+BlS4Zg306QOdOsGwYTUsoC0iTZH+2puhjh0bsFNZGZx2GrRpAyNGaJEikWZGZxZS\nt/LyMKZ29epwkaNVq7gjEpEMU7KQ2rnD738P33wTZhRs0ybuiEQkBkoWTdjy5Y2swB2uuCKsQ/HS\nS7DZZimJS0Ryj5JFE7V8Ocyalc8++xSxzTYbv7fRetg1cYdrrw136r3yCmy5ZXoCFZGcoLmhmqiz\nzoK2beGBBxpYwY03wqhRUFwMP/pRKkMTkQxpSmtwSxr8858wZQq8914DK7jjjrB+6uuvK1GICKBk\n0eR88UW4zDBuXANuuAO4/374299CoujUKeXxiUhuUrJoQsrKwpx+11wT5nyqtxEjwmIWxcWw886p\nDk9EcpiSRRPyl7+E69BXXNGAnZ98EgYNgtdeg913T3lsIpLblCyaiIkTwwwcU6c2YBaOMWOgf//Q\nd1WvqWdFpLlQsmgCli4N3U8PPgjbb1/PncePh/PPhxdegP33T0t8IpL7Ypsbysx+bWYfmlm5mR0U\nVxy5zh0uuihMBHvCCfXceeJEOOMMeOYZOOSQtMQnIk1DnGcWM4BTgIbeCSCEa9IffwyPPlrPHd95\nJ8wg+/jj8NOfpiM0EWlCYksW7v4pgJkWRGioWbNg4EB49dV6zsTx/vtw4onw8MPwi1+kLT4RaTo0\nRXmOWrsWTj8dBg+Gffetx44zZ8Jxx8G99zag30pEmqu0nlmY2Xhgu8QiwIHr3H1MfeoqKipa/7yw\nsJDCwsIURJi7Bg+G7baDSy+teZsh/fpRWlKyoWD1apg+nfwePRjUt2/6gxSRjCouLqa4uDgtdcc+\nN5SZvQZc6e5Ta9lGc0MlePXVMPpp+nTYdtuatysqLKRowoSq5b16UZSmXygRyR6pnBsqW7qhdN0i\nSYsXh3WIHnmk9kQhIpJKcQ6dPdnM5gKHAf8xs//GFUuucIcLL4Tf/AaOPTbuaESkOYlzNNRoYHRc\n7eeiBx+EL7+EJ56IOxIRaW50B3eOmDkzrEU0caJWNhWRzFOyyAFr1oRhsn/9K+y9d/L75e+yC0Vv\nvBGmoG3bdkN5QUEaohSRpiz20VDJaO6joa68MnQ/Pf001OsWxqFDYdo09VuJNFNaKa8ZGTcuzB4+\nfXo9E8Xq1XDnnaECEZFGUrLIYt9+C+edB489Bh071nPn4cPh0EPreXu3iEj11A2VpdzD9E3dusGQ\nIfXcee1a6No19FtpNlmRZkvdUM3A3/8OCxeG2cPr7bHHwiJGShQikiI6s8hCM2bAUUfBm2/CHnvU\nc+eysjBk6qGH4Mgj0xKfiOSGpjjdh0RWrw7DZG+9tQGJAsLV8E6dlChEJKV0ZpFl+vcPF7ZHjarn\n6CeAigrYbz+47bYwDbmINGu6ZtFE/ec/8PzzDRgmW+n55yE/XxNHiUjKKVlkifnzwySBTz4J7ds3\noAJ3uPFGuO66BmYaEZGa6ZpFFqiogHPPDcmiZ88GVjJuHJSWwkknpTI0ERFAySIr3HUX/PAD/PnP\njajkppvgmmsgT/+kIpJ66oaK2bRpcPPN8M470LKh/xqvvw7ffAO//W1KYxMRqaT/hsZo5cowTPbu\nu2HXXRtR0U03waBBjcg2IiK109DZGF10EaxaFW64brApU+BXv4LPP4fWrVMWm4jkvlQOnVWyyKB+\n/YZQUlIKwHffhe/3gw+GffbJZ9iwQQ2r9JRTwu3e/funMFIRaQoycp+FmS2vKw5gvrtrJZ0klZSU\nMmFC0UZlb7wBLVsWVbd53WbMgLffhscfb2xoIiK1qq2T+3N3P7C2nc1sWorjkfq4+Wa44grYbLO4\nIxGRJq62ZNEnif2T2UYiKe1J++wzGD8e7r8/hZWKiFSvxtFQ7v6FmbUws9dq2yY9YTVNX6TyaA0Z\nApdeCu3apbBSEZHq1TrW0t3LzazCzLZy92WZCqopuusuWLw4RZXNmQPPPguzZqWoQhGR2iUzMH8F\nMMPMxgMrKwvd/fLGNGxmQ4ETgTXA58B57l7XRfWcNGoU3H47nHxyPt98U1Tl/YKC/PpVeOutcP75\n0KFDagIUEalDnUNnzeyc6srdfUSjGjb7BfCqu1eY2ZBQpV9Tw7Y5O3T25ZfDjXevvJKi5bAXLgyL\nG338cVi3QkSkBpmeonwJ8IK7V6SiwUru/nLCy7dpghfLp04NieKpp1KUKADuuAPOOEOJQkQyKpnp\nPn4LzDKzoWa2V5ri+B3w3zTVHYvPP4cTTgiDlVK2aN2SJTB8OAwcmKIKRUSSU+eZhbufaWbtgNOA\nR83MgUeAJ9z9h9r2ja5zbJdYBDhwnbuPiba5Dljn7rXeWVZUVLT+eWFhIYWFhXWFHpuFC8P6Q3/+\nc5iJI2XuuSfcsd25cworFZGmori4mOLi4rTUnfR0H2bWETgLuAKYCXQF7nH3vzW4cbNzgQuBo9x9\nTS3b5cw1ix9+gMLCcFZxww0prHj5cth9d3jrLejaNYUVi0hTldG5ocysN3AeITmMBEa4+7dm1hb4\n2N13aVDDZscBtwNHunutg0pzJVmsXRuSxK67hu6nlC5Yd8st8P77mtpDRJKW6WQxAnjI3V+v5r2f\nu/srDWrYbBbQGqhMFG+7+yU1bJv1yaKiAs46K0w7/tRTKZ4tfPXqkIFefhm6dUthxSLSlGVqIsGX\ngLHAze7+SXXbNDRRRPvu0dB9s9HAgeFeuXHj0rCsxPDhcPjhShQiEpsazyzMrBNwXPQoACYTksfL\n7r6y2p3SJNvPLG67DR59FCZOhPbtU1z52rXhGsUzz4T5zEVEkpTx9SzMLA84FPgl8HNgNTDO3Yem\nIogk2s/aZPHYY3D99WGq8Z12SkMDw4fDv/8NL72UhspFpCmLffEjM/sRcKy7/zMVQSTRXlYmi7Fj\n4Zxz4LXXYJ990tBAWRnstRc88gj07JmGBkSkKcvoHdxmtitwOdAlcXt3752KAHLVlCnhgvZzz6Up\nUQA8+STssIMShYjELplLsaOBh4DngZRO+ZGrZs2C3r3hoYegR480NVJRAX/9a5iBUEQkZskki1J3\nvyftkeSIBQvC3dl/+UtIGGnz3HOQnw/HHJPGRkREkpPMfRanA3sA4wjTiQPg7lPTG9pGMWTFNYvl\ny6FXL+jTJ1zUTht3OOSQ0MjJJ6exIRFpyjI96+y+hGk+jmJDN5RHr5uNNWvCtExHHAHXXZfmxsaN\ng9LSNJ+6iIgkL5kzi8+Afdx9bWZCqjaGWM8sKirgtNPCz1GjoEWLNDfYsydcfHGY31xEpIEyfWbx\nIbA18G0qGsw17nDFFWEm2bFjM5AoXn89XBj5zW/S3JCISPKSSRZbA5+Y2RQ2vmbRLPpIhgyBCRPC\nd3h+PVc/bZCbboJBg9IwZ4iISMMl8400OO1RZIl+/YZQUlK6/vX8+TB7NvTpk89WWw1KfwBTpsDM\nmeEGDhGRLJLM4kcTMhFINigpKWXChKIq5fPmVS1Li5tuCjMStm6dmfZERJJU47KqZvafunZOZhtJ\n0owZMHkyXHBB3JGIiFRR25nFT83s+VreNyBdE100PzffDH/4A2y2WdyRiIhUUVuyOCmJ/WMbTtuk\nzJoF48fDAw/EHYmISLVqTBbN6VpFpVWrYmp4yBC47DLYcsuYAhARqZ3GZyaYPz+f3XYrYuedNy4v\nKEjjmNk5c2D06HB2ISKSpZQsIi++CG3bDmLGjAwPRrr11nBRu0OHDDYqIlI/DVr8KNPSPd3H2rWw\n775w551w/PFpa6aqBQvCYhgzZ8J222WwYRFpDlI53UdtQ2f3MLNHzewOM9vJzP5rZivN7H0zOyQV\njWeLv/0tLHOd0UQBcMcdcOaZShQikvVqPLMws0nASKAd8AfgCmAM0BO40d0PzViQaTyzWLgQunWD\nSZNgzz3T0kT1liyBPfaA6dOpcpFERCQFMrIGt5lNd/cDouefuXvX6t7LhHQmi/PPh/bt4bbb0lJ9\nzYqK4OuvYfjwDDcsIs1FpmadTVxCdXkt7zWImf0v4V6OCmAhcK67L2hsvfXx7rvhwvYnn2SmvSH9\n+lFaUgJlZeFu7YMOgsJC8gsKGDRsWGaCEBFpgNqSxV5m9gHhTu3do+dEr3dLQdtD3f3PAGbWnzBh\n4cUpqDcp7jBgANx4I2y1VWbaLC0poWhCwu0r77wDQFFmmhcRabDaksXe6WzY3VckvNycFJyt1Mfj\nj4fV7847L5Otiojkptru4J6d7sbN7EbgbGAp8LN0t1dpxQq4+mp48knIq3E8mIiIVKoxWZjZD4S1\nttcXRa8NcHdvV1flZjYeSBwXWlnHde4+xt2vB643s6uB/tTSI1NUtOGtwsJCCgsL62q+RkOGQK9e\n0KNHg6tomKVLM9ygiDQnxcXFFBcXp6Xu2kZDjQY6Ac8Ao9x9TloiCG3tDLzo7vvW8H7KRkN9+SUc\nckgYsbrTTimpMjlr1lDUvj1Fq1dXeauoVy+K0vQPLCLNV0ZGQ7n7yWa2FfAr4EEzywf+RUgcSxrb\nsJl1dffPopcnAzMbW2cyrroqzASe0UQBcMst5G+zDUW77AK28b9dfkFBhoMREamfpKb7MLM84FTg\nHuCv7n5Hoxs2ewooIFzYng383t3n17BtSs4sXn013Ffx8ccZXjbi00/hpz+FqVN1A56IZExGbsqL\nGuoBnEa4a3sS8C93n5iKhusjFcmirAwOPBBuuAF+9asUBZaMigr42c+gTx+4/PIMNiwizV1GuqHM\n7CvCKKVRQD+gLCo/CMDdp6YigEx54AHYZhs45ZQMN/zII7B6NVx6aYYbFhFJndoucBezYTRU5Sio\nSu7uR6U3tI1iadSZxeLFsPfe8MorYXbZjFm4MDQ4fjzsv38GGxYRyWA3VLZobLK47LLw8957UxRQ\nsk4/PVyjuOWWDDcsIpK5bqhDgLmV8zWZ2dlAH8LF6KJUjIjKhBkzws13MzMy1irB2LFh/idNFCgi\nTUBt9y8/AKwFMLMjgSGEKcuXATkx613l/E+DB0PHjhlseOVKuPhiuO8+aNs2gw2LiKRHbXNDtUg4\ne/gtMMzdnwaeNrPp6Q+t8Z59Fr77Di66KMMNFxXBEUfAMcdkuGERkfSoNVmYWUt3LwN+ThgRlcx+\nWWH1arjySnjoIWiZyWinTYORI0P/l4hIE1Hb1+gTwAQzWwSsBiZCuPOa0BWV1e64IywXcVTGxmwB\n5eXQr1+YfGrbbTPYsIhIetV1U95hwPbAOHdfGZUVAFtk8j6L+o6G+vprOOAAmDIFdt01jYFt6u67\nYfTocKu4pWQAgohIg2nobB3OPBN22SUsbJQxc+aEU5k33wTN9SQiWUDJohZvvgm//W1YKnXzzdMc\nWCV36N0bDj0Urr8+Q42KiNQuU2tw55yKijD90i23ZDBRADz9NHzxRfgpItIE5fyZRb9+QygpKQVg\n/vzwOPBA2HPPfIYNG5T+4JYuhR//ONz5d8QR6W9PRCRJOrNIUFJSyoQJRRuVvf46mBVVu33KXXNN\n6IJSohCRJiznk0Ws3ngDnn8ePvoo7khERNKqtuk+pDZr14Z7Ku6+G7beOu5oRETSSsmioYYOhd12\nC4saiYg0ceqGaoiSErjrrrBMqm6+E5FmIOeTRUFBPlBUQ3kauIeZCa+/Hjp3Tk8bIiJZJueHzmbc\nI4/AP/4Bb78NLVrEHY2ISI10B3dcvv02LJM6dmy4mUNEJIspWcTlzDNh++3h1lvjjkREpE66KS8O\nL70U7qv48MO4IxERybjYh86a2ZVmVmFmHeKOpUarVm1YJjWjk06JiGSHWJOFme0EHA3MjjOOOt1w\nAxx2GBx3XNyRiIjEIu5uqDuBgcDzMcdRs/ffDyOgtEyqiDRjsZ1ZmFlvYK67Z++3cHk5XHgh3Hwz\nbLdd3NFPw7zQAAAP3UlEQVSIiMQmrWcWZjYeSPyWNcCB64FrCV1Qie/VqKioaP3zwsJCCgsLUxVm\nzf7+d2jbFn73u/S3JSLSSMXFxRQXF6el7liGzppZN+BlYBUhSewEzAO6u/u31WyfkaGzQ/r1o7Sk\nJLwoLYX33oMDDyT/gAMYNGxY2tsXEUmlnB866+4fAp0qX5vZl8BB7v59HPFUKi0poWjChI0Lp0yh\nqG3beAISEckSsQ+djTh1dEOJiEh84h4NBYC77xZ3DCIiUrNsObMQEZEspmSRKBvmnxIRyUJZ0Q2V\nLfLbtKFoiy3goIM2WtQov6AgxqhEROKnZFHJnUE//BDu1v71r+OORkQkq6gbqtLEibB4MZxyStyR\niIhkHSWLSrfcAgMHavU7EZFqaPEjgA8+CDPKfvEF5Kdp7W4RkQxL5R3cOrMAGDoUBgxQohARqYHO\nLL76Cg4+GD7/HLbaKj1tiIjEQGcWqXT77XDBBUoUIiK1aN5nFt99B3vuCR99BNtvn/r6RURilPOz\nzmaNe++Fvn2VKESyzC677MLs2dm92nI26dKlC1999VVa22i+ZxYrVsCuu8Kbb8Iee6S2bhFplOh/\nxHGHkTNqOl66ZpEKw4dDYaEShYhIEprnmcXatdC1KzzzTBgJJSJZRWcW9aMzi3R54gkoKFCiEBFJ\nUvO7wF1REW7Cu+uuuCMREckZze/M4oUXwp3av/hF3JGISDOSl5fHF198AcB5553Hn//855gjqp/m\nlyyGDIGrr95ovQoRkXSzen7nFBcXk5eXx6233pqmiOqneXVDTZoECxdCnz5xRyIi9dSv3xBKSkqr\nlBcU5DNs2KCM1dFQ9b1gP3LkSDp27MjIkSMZOHBgmqJKXvNKFpqGXCRnlZSUMmFCUTXvVFeWvjq+\n/vprBgwYwMSJE3F3TjvtNO655x4AHn74YW677TYWLlxI9+7deeCBB+jcuXPSdVdatWoVTz31FMOH\nD+fss89m6tSpHHTQQfWuJ5WaTzfUhx/Cu+/COefEHYmI5KiKigpOOOEEdt11V+bMmcO8efM49dRT\nAXjuuecYMmQIo0eP5rvvvqNnz56cdtppDWrn6aefZsstt6Rv374cc8wxjBgxIpUfo0GaT7IYOhQu\nv1zTkIs0MRMmhEuQyTwmTGhcW++88w7z589n6NCh5Ofn07p1a3r06AHAAw88wDXXXENBQQF5eXkM\nGjSI6dOnM3fu3Hq3M3LkSE499VTMjNNPP51Ro0ZRXl7euOAbKbZkYWaDzexrM5saPY5LW2OzZ4dR\nUBdfnLYmRCQevXqBe3KPXr0a19bcuXPp0qULeXlVvzpnz57NgAED6NChAx06dKBjx46YGfPmzatX\nG19//TWvvfYap59+OgC9e/dm9erVvPDCC40LvpHivmZxh7vfkf5W7oDzz4ett057UyLSdO28887M\nmTOHioqKKgmjc+fOXH/99Q3ueqo0cuRI3J0TTzxx/UXxNWvWMGLECHr37t2ouhsj7mSR/vGrixbB\nY4+FaxYikrMKCvKp7kJ0KM9MHd27d2f77bdn0KBBFBUV0aJFC9577z169OjBRRddxJ/+9Cf2339/\n9tlnH5YtW8b48eP59a9/nXR8EJJFUVERF1100fqyyZMn07dvX77//nvat29fr/pSJe5kcZmZnQW8\nC1zp7stS3sK994ahsjvskPKqRSRzUjG0tbF15OXlMWbMGPr370/nzp3Jy8vj9NNPp0ePHpx88sms\nXLmSU089lTlz5rDVVltx9NFHr08WydxnMXnyZObMmcMll1xCx44d15efeOKJ7LHHHjzxxBNccskl\njfoMDZXWiQTNbDywXWIR4MB1wNvAInd3M7sR2N7dz6+hHh88ePD614WFhRQWFtYdwMqVYRrySZPC\nXFAikhM0kWD9VB6v4uJiiouL15ffcMMNKZtIMCtmnTWzLsAYd9+vhvcbNuvsPfeE4Q9PP93ICEUk\nk5Qs6icTs87G1g1lZp3cfUH08ldAai8qrFsX1tf+979TWq2ISHMU5zWLoWZ2AFABfAVcVPvm9TRq\nFOy+O3TvntJqRUSao6zohqpLvbuhKipgv/3CmcWxx6YvMBFJC3VD1Y8WP2qoF1+EVq3gmGPijkRE\npElomsnills0DbmISAo1vWTxxhswbx7U80YYERGpWdNLFpXTkLeM+35DEZGmo2kli48+gnfegXPP\njTsSEZGNaFnVbDJ0KPTvD5ttFnckIiIbSXZZ1REjRtCyZUvatWtHu3bt6Nq1K/fff3+ao6tbzvfV\nDOnXj9KSEigtDYsbHXoojB9PfkEBg4YNizs8EUmR9X/rm6jP33oq6mio+gwF7tGjB6+//joA06dP\n58gjj+Twww9n//33T1d4dcr5ZFFaUkJR4oomb74J1GeRRBHJBVX+1iNFGa4jE8uqJjrggAPYe++9\nmTlzZqzJoml1Q4mIpFGmllVNNGXKFGbNmsXBBx/c6LoaQ8lCRHJbBtdVzdSyqm+99RYdOnSgXbt2\nHHbYYZx11ll07dq1UbE3lpKFiOS2DK6rmollVQEOP/xwlixZwvLly1mwYAEffvgh1113XaNibywl\nCxGRJCUuq7qpzp0788ADD7BkyRKWLFnC999/z4oVKzjssMMa1eY222xDnz59GDNmTKPqaaycv8Cd\nX1BQ7cWpfC12JNKkpOJvvbF1ZGJZVdh45NTixYt59tln6datW73rSaWcTxYaHivSPKTib72xdaR7\nWdVKb7/9Nu3atQOgbdu2/OIXv+Cuu+5qVOyN1TSnKBeRnKYpyutHU5SLiEhWULIQEZE6KVmIiEid\nlCxERKROShYiIlInJQsREalTzt9nISJNT5cuXep1X0Jz16VLl7S3Eet9FmbWH7gEKANecPdBNWyn\n+yxEROqpSdxnYWaFwInAvu6+L3BbXLGkW3FxcdwhNEoux5/LsYPij1uux59KcV6zuBgY4u5lAO6+\nKMZY0irXf+FyOf5cjh0Uf9xyPf5UijNZFABHmtnbZvaamcW7soeIiNQorRe4zWw8sF1iEeDA9VHb\n7d39MDM7BHgS2C2d8YiISMPEdoHbzF4EbnH3CdHrz4BD3X1xNdvq6raISAOk6gJ3nENnRwNHARPM\nrABoVV2igNR9WBERaZg4k8UjwMNmNgNYA5wdYywiIlKLnFjPQkRE4pXV032Y2XFm9omZlZjZ1XHH\nUxMz+8rM3jezaWb2TlTW3szGmdmnZvaSmW2VsP01ZjbLzGaa2TExxPuQmS00sw8Syuodr5kdZGYf\nRP8+GVvGq4b4B5vZ12Y2NXocl43xm9lOZvaqmX1kZjPM7PKoPCeOfzXx94/Kc+X4tzGzydHf6gwz\nGxyV58rxryn+9B9/d8/KByGRfQZ0AVoB04G94o6rhli/IIzsSiy7Bfh/0fOrCfeUAOwDTCN0Ae4S\nfUbLcLw/BQ4APmhMvMBk4JDo+YvAsTHGPxj4YzXb7p1N8QOdgAOi51sAnwJ75crxryX+nDj+UVtt\no58tgLeB7rly/GuJP+3HP5vPLLoDs9x9truvA0YBJ8UcU02MqmdpJwEjoucjgJOj572BUe5e5u5f\nAbMInzVj3H0S8P0mxfWK18w6AVu6+5Rou5EJ+6RVDfFD+HfY1ElkUfzuvsDdp0fPVwAzgZ3IkeNf\nQ/w7Rm9n/fEHcPdV0dM2hC9RJ0eOP9QYP6T5+GdzstgRmJvw+ms2/FJmGwfGm9kUM7sgKtvO3RdC\n+AMDto3KN/1c88iOz7VtPePdkfBvUikb/n0uM7PpZjY8oRsha+M3s10IZ0hvU//fl2yKf3JUlBPH\n38zyzGwasAAYH31h5szxryF+SPPxz+ZkkUuOcPeDgOOBS82sJxuyfaVcG0mQa/H+A9jN3Q8g/BHd\nHnM8tTKzLYCngAHR/9Bz6velmvhz5vi7e4W7H0g4o+tuZj8mh45/NfHvQwaOfzYni3lA54TXO0Vl\nWcfd50c/vyPcP9IdWGhm2wFEp3zfRpvPA3ZO2D1bPld9482qz+Hu33nU+Qo8yIauvayL38xaEr5o\nH3P356LinDn+1cWfS8e/krsvB4qB48ih418pMf5MHP9sThZTgK5m1sXMWgOnAs/HHFMVZtY2+l8W\nZrY5cAwwgxDrudFm5wCVXwrPA6eaWWsz2xXoCryT0aADY+M+znrFG52qLzOz7mZmhPtkniNzNoo/\n+gOv9Cvgw+h5Nsb/MPCxu9+dUJZLx79K/Lly/M3sR5VdNGa2GXA04bpLThz/GuL/JCPHPxNX7xtx\n1f84wmiLWcCguOOpIcZdCSO1phGSxKCovAPwchT/OGDrhH2uIYxKmAkcE0PMjwPfEG6GnAOcB7Sv\nb7zAT6LPPAu4O+b4RwIfRP8Wowl90FkXP3AEUJ7wOzM1+j2v9+9LlsWfK8d/3yjm6VG810XluXL8\na4o/7cdfN+WJiEidsrkbSkREsoSShYiI1EnJQkRE6qRkISIidVKyEBGROilZiIhInZQsRESkTkoW\nktPMrDyav39a9PP/ReWvmdlBKW5rgJnlJ7z+j5m1S0G9j5jZF2bWr577vWpmP6T6c4pUJ85lVUVS\nYaWHSRwz4QrgMaAUwN1PSGHdV7n7M/XZwd2PMrNXUxiDSI10ZiG5rro5/DfewOxoM3vTzN41s39F\n83kda2ZPJmzTy8yej57/w8ze2WQlsv7ADsBrZvZKVPalmXWInv8x2v4DMxsQlXUxs4/NbJiZfWhm\nY82sTRLxPhLF8JaZfRbF9lBU18P1/fwiqaBkIblus026ofomvmlmHYHrgZ+7+8HAe8AfCfMAdY8m\nYwP4LWGBLYBr3b07sD9QaGbd3P1vhFk5C93959F2HrVxEGHyuUOAw4ELzWz/aJuuwN/cvRuwDOiT\n5Ofa2t0Pj2J9Hrjd3fcB9jOz/ZKsQyRl1A0luW5VHd1QhxGWxnwjml2zFfCmu5eb2VjgRDN7Gvgf\nYGC0z6lmdiHh76NTtP+HVJ2pt9JPgWfdvRTAzJ4BegJjgC/dfUa03XuEpS2TMSb6OQNY4O4fR68/\niur4oLqdRNJFyUKaOgPGufsZ1bz3L+AywhKtU9x9pYXV364EfuLuy83sESC/mn2TtSbheXk96qrc\nr2KTOirQ363EQN1Qkuvq6rN/GzjCzHaH9euP7BG9NwE4CLiQDV1Q7YAVwA/RYji/TKhrefT+pm1P\nBE42s/xoTZNTorJk4kuGrktI7PQ/FMl1+WY2lfCF6sBYd782eo67LzKzc4EnoovLTriGMcvdK8zs\nP4TrDWdH239gZtMJc//PBSYltPUgMNbM5kXXLSrbmGZmjxIW7HJgmLu/b2ZdaNjynLUt8ak1BSQW\nWs9CJGZRV9d/3P3pBuz7GnClu09NfWQiG6gbSiR+y4D/bchNeYSVGtelJSqRBDqzEBGROunMQkRE\n6qRkISIidVKyEBGROilZiIhInZQsRESkTv8fo53TTSiF4XAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a2c6250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(elev_mec,smb['cellA'],'bs-',label='cell A')\n",
    "plt.plot(elev_mec,smb['cellB'],'rs-',label='cell B')\n",
    "plt.xlabel('Elevation [m]')\n",
    "plt.title('Surface mass balance')\n",
    "plt.ylabel('SMB [m/yr]')\n",
    "plt.legend(loc='lower right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Downscaling options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose we would like to downscale SMB to a CISM grid cell *P* which is horizontally located exactly midway between cell A and cell B, with some arbitrary elevation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "elev_P = 500. #m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In reality we would consider bilinear interpolation but in this case we can simplify to a linear interpolation, I think without loss of generality. There are several options how to do the downscaling. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Method 1: horizontal then vertical\n",
    "\n",
    "Horizontally interpolate SMB at each elevation class, then vertically interpolate to the target elevation. \n",
    "\n",
    "Number of linear interpolations required: 2x10 + 1 = 21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.5628565450650003"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smb_P = {}\n",
    "from scipy.interpolate import interp1d\n",
    "smb_interpolated = 0.5*(smb['cellA']+smb['cellB'])\n",
    "y_interp = interp1d(elev_mec, smb_interpolated)\n",
    "smb_P['method1'] = float(y_interp(elev_P))\n",
    "smb_P['method1']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Method 2: vertical then horizontal\n",
    "\n",
    "First, vertically interpolate SMB to the target elevation, then do the horizontal interpolation. \n",
    "\n",
    "Number of linear interpolations required: 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.5628565450649998"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_interpA = interp1d(elev_mec, smb['cellA'])\n",
    "yA = y_interpA(elev_P)\n",
    "y_interpB = interp1d(elev_mec, smb['cellB'])\n",
    "yB = y_interpB(elev_P)\n",
    "smb_P['method2'] = 0.5*(yA+yB)\n",
    "smb_P['method2']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same answer as with option 1, as we would expect since these are linear operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Method 3: regression on AB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider accumulation and ablation separately. In each CLM grid cell, derive regression parameters a and b for a least-squares regression of ablation with height (this is the method that Noël2016 uses:  http://www.the-cryosphere.net/10/2361/2016/). Then interpolate regression parameters a and b to the target location, call these a_P and b_P. Compute now downscaled ablation as\n",
    "```\n",
    "AB_P = a_P + b_P * elev_P\n",
    "```\n",
    "For accumulation no elevation correction is required, as this field does not always scale the same way with height (?)\n",
    "\n",
    "In the test data, ablation (AB) shows a strong correlation with elevation up to about 1200m, after which it levels out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10fb17b50>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEZCAYAAABmTgnDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOX5///XFRaDCi6IK4JajGsRUaNiLakVd1x+Sota\nF/qr2KoVP1orVq1xq7jVpdUKKgqtdUNFcUFxCe64IO4WqJZFES0gIDQIyfX94z6BkEwyM8nMnJnJ\n+/l4zGNmzpxz31cGyMV97s3cHRERkeaUxB2AiIjkPyULERFJSslCRESSUrIQEZGklCxERCQpJQsR\nEUlKyULaNDO728wub+bzWjPbroVln2BmE1senUj+ULKQNsPMqsxsoZl1SOOylCYimVnPKLGs/jfl\n7v9090PSDlQkDylZSJtgZj2BHwG1wJHpXJrGeZ7G+SIFRclC2oqTgdeBe4BTG3zWzcyeNbMlZvai\nmfVIVICZHWZmU81ssZnNMrNL6308OXr+NipnbzM7xcxernd9PzN708wWmdkUM9u33mcvmtnlZvZK\ndP1EM9s4Ez+4SCYoWUhbcTLwD+CfwMFm1q3eZycAlwFdgfeAe5so4zvgJHffADgc+LWZ1bVSfhw9\nd3H3Lu4+JXrvAGa2EfAEcFNUz43Ak9HxOscDpwDdgHWA37XwZxXJOCULKXpm9iOgB/Cgu08FZhIS\nRJ0n3f1Vd18JXATsa2ZbNSzH3V9y94+i1x8C9wP9G1bXRBiHA9Ojfoxad78f+BQYWO+cu9393+6+\nAngQ6JP2DyuSJUoW0hacDDzr7oui9/cR/gdfZ07dC3dfBiwEtmxYSHRr6QUz+9rMvgVOBzZJMYYt\ngVkNjs0C6ielr+q9Xg6sn2LZIlnXPu4ARLLJzEqBnwElZjYvOrwOsIGZ9Y7eb13v/PWBjYEvEhR3\nL3ALcLC7rzSzGwm3lCD5qKkvgWMbHOsBPJ3qzyISJ7UspNgdA6wCdgJ2ix47Ai8TWhwAh0Wdzx2B\nK4DX3f3LBGWtDyyKEkU5a9/K+oYw0uoHTcTxFLC9mQ02s3Zm9vMopgmt+/FEckPJQordycBod//C\n3b+uewC3En7ZtyN0elcCC4DdgV/Uu75+i+EM4AozWwxcDDyw+iT3/wFXAa9GcznK6wfh7guBIwid\n1v+Nng+vd2tMG8tIXrO4Nz8ysw2AO4FdCf8z+2W9kSQiIpIH8qHP4mbgKXcfZGbtgXXjDkhERNYW\na8vCzLoA77p7U/d5RUQkD8TdZ7Et8N9oMbepZjbKzDrFHJOIiDQQd7JoD/QFbnX3voSx5cPjDUlE\nRBqKu89iLjDH3d+O3o8DLmh4kplppIiISAu4e0YWt4y1ZeHu84E5ZlYWHfop8HET5xbs49JLL409\nhrYafyHHrvjjfxR6/JkUd8sC4Gzg3miPgc+AITHHIyIiDcSeLNz9PWCvuOMQEZGmxd3B3SZUVFTE\nHUKrFHL8hRw7KP64FXr8mRT7DO5UmJkXQpwiIvnEzPBi6OAWEZHCoGQhIiJJKVmIiEhSsY+Gaq0R\nQ4dSPX16o+OlZWUMHzUqhohERIpPwSeL6unTqZw8udHxytyHIiJStHQbSkREklKyEBGRpJQsREQk\nKSULERFJquA7uEvLytZ0Zi9ZAh9/DOXllJaVNXOViIiko/iW+xgwAE44AYZo8VoRadsyudxH8SWL\nyZPhV7+CTz+Fdu2yG5iISB7T2lDN+fGPYfPN4cEH445ERKRoFF/LAuCZZ+C88+D996Gk+PKhiEgq\n1LJI5qCDoFMneOyxuCMRESkKxZkszODii+HKK6EAWk4iIvmuOJMFwMCBsHIlTJwYdyQiIgWveJNF\nSQlcdBFccYVaFyIirVS8yQLguONgwQKoqoo7EhGRglbcyaJdO7jwwtB3ISIiLVbcyQLgxBPh3/+G\n11+POxIRkYIV+zwLM/sPsBioBVa6e3mCc9KbZ9HQ7bfDE0+Eh4hIG1FUy32Y2WfAHu6+qJlzWpcs\nqqvhBz8IyWL33VtejohIASm2SXlGtuMoLYXzz4errspqNSIixSpfWhbfAjXAKHe/I8E5rWtZACxb\nFloXzz8Pu+zSurJERApAJlsW+bCfxX7uPs/MugGTzOwTd3+l4UmVlZWrX1dUVFBRUZFeLeutB+ec\nA1dfDf/4R+siFhHJQ1VVVVRlaapA7C2L+szsUmCpu/+5wfHWtywgbI603XbwxhvQq1fryxMRyWNF\n02dhZuua2frR6/WAg4APs1Zhly5w5pkwYkTWqhARKUaxtizMbFvgUcAJt8TudfdGv8kz1rKAMKN7\n++1h2jTo0SMzZYqI5KGiGjqbiowmC4ALLggd3n/9a+bKFBHJM0oWrTV/Puy0E3z0EWyxRebKFRHJ\nI0XTZxGbzTaDk06CG26IOxIRkYLQNlsWAHPnQu/eMH06bLJJZssWEckDallkQvfuMGgQ3Hxz3JGI\niOS9ttuyAPjsMygvh5kzYcMNM1++iEiM1LLIlO22g8MOg1tvjTsSEZG81rZbFgCffAL9+4dWxvrr\nZ6cOEZEYqGWRSTvtBBUVMHJk3JGIiOQttSwA3nsPDj007KjXqVP26hERySG1LDJtt91gzz1h9Oi4\nIxERyUtqWdSZMiUMpZ05Ezp2zG5dIiI5oJZFNuy9N+y4I/z973FHIiKSd9SyqGfEwIFUP/dcmHth\na5JxaVkZw0eNynr9IiKZVGw75eWN6qVLqayuhpdeWut4ZTzhiIjkDd2GEhGRpJQsREQkKSULERFJ\nSslCRESSUgd3PaVlZWs6s91h6lTYYgtKy8pijEpEJH4aOtucDz+En/wE3nkHevTIff0iIq2gSXm5\nsuuucM45MHRoaGmIiLRRShbJ/P738PXXcM89cUciIhIb3YZKxXvvwYEHwrRpsNVW8cUhIpKGorsN\nZWYlZjbVzB6PO5aEdtsNzjoLTj9dt6NEpE3Ki2QBDAM+jjuIZl14IcyZA//4R9yRiIjkXOzJwsy6\nA4cBd8YdS7M6doS774bzzoN58+KORkQkp2JPFsCNwPlA/t/f6ds3jIz6zW90O0pE2pRYJ+WZ2eHA\nfHefZmYVQJMdMZWVlatfV1RUUFFRke3wErvkEthjD7j/fjj++HhiEBFJoKqqiqqqqqyUHetoKDP7\nE/ALYBXQCegMPOLuJzc4L97RUA299RYccQS8/z5stlnc0YiIJJTJ0VB5M3TWzPoD57n7kQk+y69k\nATB8eNiCddy4uCMREUlImx/VM3ToCKZPr250vKyslFGjhmev4spK2H13eOihsHe3iEgRy5tk4e6T\ngcnpXjd9ejWTJ1cm+CTRsQwqLQ2jo445BioqoFu37NYnIhKjfBgNVbj22QdOPBHOPjvuSEREskrJ\norWuuCIsZf7oo3FHIiKSNUoWrdWpE4weDWeeCQsWxB2NiEhWKFlkwn77wc9+FpYzFxEpQnnTwd1S\nZWWl1HVmf/ddWCB2r73qjufQVVeFBQcnTICBA3Nbt4hIluXNPIvmpDPP4txzYdGiMFAp5yZPDh3e\nH3wAG20UQwAiImsU5aS85qSTLJYuhZ12gvvug/33z3JgiZx1FixbFlO2EhFZo+j2s8ikzp3hxhvh\njDNg5coYAhgxAqqq4OmnY6hcRCQ7ii5ZABx3XNjQ7uabY6h8/fXhrrvCRkmLF8cQgIhI5hXdbag6\nM2eGOXPvvgtbb52lwJrz619DTQ3ccUcMlYuIqM8iZZddFhaGffjhLASVzJIl8MMfhmRx0EExBCAi\nbZ36LFJ0wQUhWTz1VAyVd+kSEsXQoaHXXUSkgBV1ywLg2WfDHaGPPgqTrXPuV7+CDh3gb3+LoXIR\nact0GypNP/sZ7LBDWMYp5xYvhl13hTFj4IADYghARNoqJYs0ffFFmFz96qshaeTcU0+F+Rfvvx9G\nS4mI5ICSRQvceCM8+SRMmgSWka8uTaeeGiaB/OUvMVQuIm2RkkULrFoFe+wBF14IgwdnKLB0LFrE\niO7dqd5+e9hww7U+Ki0rY/ioUTEEJSLFTNuqtkD79qGPedAgOPRQ2GCDHAew0UZUb7stle+91+ij\nyhyHIiKSrqIeOttQv35w2GFwySUxBbDJJjFVLCLSOm0qWUBYuumBB8LmdiIikpo2lyy6doWrr4bf\n/CasxiEiIsm1uWQBYWBS+/Zw551xRyIiUhia7OA2s1tSuH6Ju1+cwXhyoqQkdHYfeCAccwxsumlu\n6i0tK1vTmT1nTlg/apddKC0ry00AIiIt1OTQWTObBfwxyfXD3X2nFldutg7wEtCRkLjGuftlCc5r\n9dDZRM47DxYsgHvuyXjRyS1bBr16wTPPQO/eMQQgIsUuJ/MszOwcd78pSSBJz0kagNm67r7czNoB\nrwJnu/ubDc7JSrJYuhR23hnuvRd+/OOMF5/cn/8Mr70G48bFULmIFLucrDrr7jeZWTsz+7/mzmlt\nAO6+PHq5DqF1kbNZgnW76v3mNzHtqvfrX4c1SBLMvRARySfNdnC7ew1wQjYDMLMSM3sX+AqY5O5v\nZbO+ho49Fnr0CEkj59ZdF37/e7j88hgqFxFJXdLlPszsRqAD8ACwrO64u2d0poKZdQHGA2e5+8cN\nPvNLL7109fuKigoqKioyVvfMmbD33mFXvR49MlZsapYvD30XTz0FffrkuHIRKSZVVVVUVVWtfn/Z\nZZflbm0oM3sxwWF394yvt21mlwDL3P3PDY5npc+ivssug2nT4NFHs1pNYjfdBC+9BI88EkPlIlKs\ncrqQoJm1i25HZZyZbQKsdPfFZtYJeAYY4e5PNTgv68miujrsgnrjjXDEEVmtqrH//Q9+8AO1LkQk\no3KdLD4DHgZGu/snmai0Xtk/BMYQ+k5KgAfc/aoE52U9WQAcdtgIXnihmr32gnbt1hwvKytl1Kjh\n2a38pptg8uSYmjYiUoxyversbsBg4C4zKwFGA/e7+5LWVu7uHwB9W1tOpixfXs2KFZW88krDTyqz\nX/npp8O114aOk913z359IiJpSLrch7svdfc73L0fcAFwKTDPzMaYWa+sR9hWdOoEF1wQOk9ERPJM\n0mQRzbU40sweBW4CbgC2AyYATzV7saRn6FB4663QuhARySOpLCQ4AzgKuM7dd3f3P7v7fHcfB0zM\nbnhtTF3rorIy7khERNbSZLIws+PNrCvQ293/f3d/reE57n52VqPLE8uWJT8nY047Dd5+WxtuiEhe\naa6DuwfwENDBzJ4HngbezMmwpJiUlZXSsDN7/nz4+ONSZsyA7bfPQRCdOsHw4aHv4rHHclChiEhy\nqQyd7QwcCBwClAOfEG4/PePu87MeIbkbOtuUO+8MK3JMngzbbpuDCqurw6zuxx6DPfbIQYUiUoxy\nOs8iQeU7A4cCB7n7wZkIIoU6Y2/Q3Hor3HBDmGjdvXsOKvzrX+HZZ+Hxx3NQmYgUo5wnCzPrDWxD\nvdtW7p6ztSnyIVkAXH89jBoVWhhbbJHlyupaF+PHw557ZrkyESlGuZ7BPRroDXwE1EaH3d1/mYkA\nUpEvyQLgyivhvvugqgq6dctyZbfeChMnwoQJWa5IRIpRrpPFx+6+cyYqa6l8ShYAF10UlnF64QXY\naKMsVqTWhYi0Qk42P6rn9aifQiJXXgkHHAAHHxy20c6a0lK48ELNuxCR2KXSsugPPE7YnGgFYITb\nUDnbODrfWhYA7nDWWWGTu4kTYf31s1TRihWhdfHII7DXXlmqRESKUa5vQ80EzgU+YE2fBe4+KxMB\npCIfkwVAbW2YQ/f55/Dkk2GKRFbcdluo4Mkns1SBiBSjXCeL191930xU1lL5miwAamrg5JNhwYIw\nLWKddbJQyYoVYUbguHFQXp6FCkSkGOU6WdwGbEhYOHBF3fG2OHS2KatWweDB4fmhh6BDhyxU8re/\nwRNPqHUhIinLdbK4O8HhNjt0tinffw/HHhtuRf3zn9A+lZ1C0lHXunjoobBhuIhIErHO4I5DISQL\nCCNdjzwSNt8c7rkHSlIZa5aO228PM7qf0srwIpJcTobOmtnQFAJJek5bUloapkTMng2//nUYMZVR\nQ4bARx/BG29kuGARkeY12bKI9t7+XXPXApe7+y7ZCKxBLAXRsqizdGmYg7HnnnDzzWAZyeuRkSND\nRnr66QwWKiLFKCe3oZroq2hosbufk4lAmlNoyQJg8WL46U/D5L1rrslgwvj++9B38cADsM8+GSpU\nRIqR+iwKxIIF8JOfwDHHZHhr7VGjwiS9idqoUESaluvlPqSFunaF556DBx+Eq6/OYMGnngqffAKv\nv57BQkVEmqaWRQ58+SX07w+bbDKCddapbvR5WVkpo0YNT6/QUaPg4YfhmWcyFKWIFJtMtiwyPRsg\nLWbWHRgLbEZYSuQOd78lzpiyYcst4fnnoaysmhUrKhOckehYEqeeCn/6E7z2GvTr17oARUSSSHob\nysy6mtlfzGyqmb1jZjebWdcM1b8KODcaUbUvcKaZ7ZihsvNKjx6w224ZLLBjx7BWulakFZEcSKXP\n4n7ga+BY4DjgG+CBTFTu7l+5+7To9XeE/b23ykTZ+SjjCw2ecgrMmAGvvprhgkVE1pZKstjC3a9w\n98+jx5WE20YZZWbbAH2AKZkuu2jVtS4yOtRKRKSxVPosnjWzwcCD0fvjgIz2qprZ+sA4YFjUwmik\nst7tloqKCioqKjIZQuE65RS46qrQuthvv7ijEZEYVVVVUVVVlZWym5uUtxRwwkzt9Vizl0UJ8J27\nd8lIAGbtgSeAp9395ibOKejRUHWGDh3B9OlhNJQ7fPABrLceHHpoC0ZD1XfXXXD//TBpUoYiFZFi\nUFST8sxsLPBfdz+3mXOKIlk0tGAB9O0Lt9wCRx3VioJWroQddoCxY+FHP8pYfCJS2HKeLMxsI2B7\noLTumLu/1OrKzfYDXiLswufR4w/uPrHBeUWZLCCsCXjkkeF5u+1aUdDo0WFt9Oeey1hsIlLYcr2f\nxa+AYUB3YBqwD/C6ux+QiQBSUczJAsJig2PGhCkTpaXJz09o5UrYccewNvr++2cyPBEpULle7mMY\nsBcwy91/AuwOfJuJyiU4++zQqjinNUsydugAF1+seRcikhWpJItqd68GMLN13P1TYIfshtW2mIW7\nSC+8APfe24qCfvEL+M9/4KVW3yEUEVlLKslirpltCIwHJpnZY8Cs7IbV9nTpAuPGhdbFxx+3sBC1\nLkQkS9IaDWVm/YENgInu/n3Wompcb1H3WdQ3ejRcfz28+Sasv34LCli1KvRd3HVXWL1QRNqsoho6\nm4q2lCwg7J76/ffwj3+0cNOke+4Jw2hfeCHToYlIAdF+FkXu1lvDhL2RI1tYwC9+ETYCnzw5o3GJ\nSNullkWemj49rN4xcSLssUcLChgzJrQwXnwx06GJSIFQy6INKCuD226DQYNg0aIWFHDiiTB3LmRp\nnRgRaVvUsshzZ58Ns2bB+PEt6L8YOzb0mCthiLRJalm0IddfD/Pnh+e0nXBC2NNVt6JEpJXUsigA\ns2dDeTk89FB6K3mMGDqU6pdegnnzoE+f1U2T0rIyho8alaVoRSRfFM0e3JKaHj3g7rvh+OPhnXdg\nsxS3nqqePp3Kf/0rvKk3q7sy8yGKSJHTbagCceihYZ+jE06Ampq4oxGRtkbJooBcdlnYNEm7qIpI\nrilZFJD27cOWFXfdBc+0ZmPbVasyFpOItA1KFgVm881DwjjlFJgzp4WFvPcefKtV5kUkdergLkD9\n+4fVaX/+8zCFomPHxOeVlpU17sx2p3TBAjj44NA82XDD7AYrIkVBQ2cLVG1t2Le7Vy+48cY0L3YP\n2eaNN5QwRIqYJuUJJSVh+afx4+Hhh9O82Axuugn23ju0MBYvzkqMIlI81LIocG+9BYcfHvbv7tUr\nzYvdYdgwmDIFnn0WNtggKzGKSDzUspDV9toLLr0UjjsO/ve/NC82g5tvDi2Mgw5SC0NEmqSWRRFw\nD5P11lsP7ryzhQWohSFSdNSykLWYwahR8MoroR+jRQWohSEizYg9WZjZXWY238zejzuWQta5M4wb\nB7/7XdhlL211CaO8XAlDRBqJPVkAdwMHxx1EMdh1V7jhhrBh0tKlLSjADG65JSQMjZISkXryos/C\nzHoCE9y9dxOfq88iDaedBkuWwP33t2DDJAh9GGefHYZaPfOM+jBECpT6LKRZt9wS9vC+9dYWFlDX\nwthrL7UwRAQooOU+KisrV7+uqKigoqIitljyXadOof9i333DHaXy8hYUUpcwzj57zdIgamGI5LWq\nqiqqsrSNsm5DFbEBA0bw8svV7LEHdOiw5nhZWSmjRg1PrRB3+O1v4e23lTBECkwx7pRn0UMyaOXK\nalasqOS11xp+Upl6IWbwl7+EhKEWhkibFXufhZn9E3gNKDOz2WY2JO6YpIG6hLHnnnDIIerDEGmD\nYk8W7n6Cu2/p7uu4ew93vzvumCSBuoSxxx5KGCJtUOzJQnJv3rzQFZE2JQyRNkvJog364gs47DCY\nO7cFF9cljL59lTBE2pB86eCWLCgrKyVRZ3avXqX06BF+3197bdiiNa3Je2bw17/CWWeFhDFxojq9\nRYpcXgydTUZDZ7Nj2jQ49VTo3j0sRLjllmkW4B4SxtSpYZRUly7ZCFNEWkgzuCUj+vSBN98MXRB9\n+sDf/55mX0ZdC6Nv3zCsdsmSrMUqIvFSy0KA0Dg49VTYdlsYORI23zyNi93hzDPh3XfVwhDJI5ls\nWShZyGorVsAVV8Add8CNN8Lxx6fRlxEljBHjxlFdVgbt1+4OKy0rY/ioUZkPWkSapGQhWfX226HT\ne4cd4G9/g802S/FCdyq32orKefMafVTZvz+VWVqzRkQSU5+FZNWee8I774Rk0bs3PPBAiheawfbb\nZzU2EYmHkoUkVFoKV18Njz8OlZVhQ6VvvknhwhZtoCEi+U7JQpq1996h33rbbUMrY9y4Fhb03/+2\ncNq4iOQDJQtJqrQ0TN575BG46CIYPDj87k/L559D//4wZUpWYpTiss0222BmeqT42GabbbL+Z6IZ\n3JKyffcNE/kuvhh++EO47TY45pi1zyktK0u4AHppr17Qrx8ce2wo6E9/Uv+GNGnWrFloUEvqLAe3\nfzUaSlrklVdgyJCwC98tt0DXrileuHw53Hwz3HBDGJt7ySWw6aZZjVUKTzSKJ+4wCkZT35dGQ0ns\nfvQjeO896NYt9GU8/niKF667Llx4IXzyCZSUwE47wZVXwrJlWY1XRFpHLQtptZdeCq2M/fYDsxHM\nmlXd6Jwmt3L9979DR8jLL4dhV0OGNJrQJ22PWhbpyUXLQslCMmLZMhg+HEaNquT77ysbfd6/fyVV\nVY2Pr/bmm/D738PXX8OIETBwoIbhtmFKFunRbSgpGOutF7a52HHHFhZQXg4vvgjXXQd/+INGTknR\nKSkp4bPPPgNgyJAh/PGPf4w5ovQoWUhGbbRR4uMLF4a+7WaZweGHh86QU08NI6d+9jOYOTPTYYrk\nXLojlqqqqigpKeG6667LUkTp0c1hyYlZs8Kgp/JyGDAgPPr2DX3cjbRrB7/8ZZjQcdNNsM8+Gjkl\nDB06gunT0+gPy1IZLZXubbWxY8fStWtXxo4dy/nnn5+lqFKnZCE5sfvuMGECVFXBpElw0klh+ZAD\nDliTPBrNK1p33XBL6rTTwoipnXeGc86B//u/cN9L2pTp06uZPLkywSeJjmWvjLlz5zJs2DBefvll\n3J3jjz+eW265BYDRo0dz/fXXM3/+fMrLyxk5ciQ9evRIuew6y5cvZ9y4cdx5552cfPLJTJ06lb59\n+6ZdTibpNpRkVFlZKf37VzZ6lJWV0rlz6Le+5ZYwcvbdd+HQQ0NXRXk5lJWFbTHGj4dvv61XaLdu\nYW7GlCnw4YfhxDvugFWrYvs5pW2qra3liCOOYNttt2X27Nl88cUXDB48GIDHHnuMESNGMH78eL75\n5hv2339/jj/++BbV8/DDD9O5c2cGDRrEQQcdxJgxYzL5Y7SIRkNJXqithfffD62OSZPg9dfDLPG6\nVsfee0OHDtHJ9UdOXXMNIx5/nOoZMxqVmcs9NEYMHUr19OmxxlBMEo3uqaiobKZVkOh4IonPTTpa\nL/LGG29w1FFHMW/ePEoa3EM97LDDGDRoEEOGDAFCYuncuTOffvopW2+9NSUlJcycOZPtttuOIUOG\nsPXWW3P55ZcnrGfAgAH07t2bG264gfvvv59hw4bx5Zdf0q5du4Tn52I0VOy3oczsEOAmQivnLne/\nJuaQJAYlJWFr1z594Pzz4X//g1dfDYnj7LPDdIz+/euSRzk7vPAi9vRTcMEFLJw+k2tXrmhU5vH/\nmp2z+N+b8Bz3ffV5rDGkqtgSW//+4fZmKioqYPLkltc1Z84cevbs2ShRQFiiZNiwYZx33nlA6KMw\nM7744gu23nrrlOuYO3cuL774IiNGjADgyCOPZOjQoTz55JMceeSRLQ++lWJNFmZWAvwV+CnwJfCW\nmT3m7p/GGZfEr1MnOPDA8LjmmtC/8fzzIXlcdx24GwMGHM5BFx7CN7/aFBIki86LFsH114cO83bt\nwmS/utf1H4mOp3luzXeN6wdYsXxVGAZWd11JSezzRwopseWbrbfemtmzZ1NbW9soYfTo0YOLL764\nxbee6owdOxZ3Z+DAgatbCytWrGDMmDEpJYvDdiyn41dftyqGROJuWZQDM9x9FoCZ3Q8cBShZyFq6\ndQuDowYPDiudT58eEsd9D7ZjcXVnYGHji1bV8M6T82jnNZT4qui5waN21erXtvpYOG7Ra/MarDZ6\nXbsqel5z3Gpr+MF38xPGvcuSOdRu3BWrrYFVqzB3fK2k0371e6s73r7B+5Y8N/PZRgsTLxncedHC\nkJlLStYktrqykr3O5LkJh8iF/rBEt5DC8dS0tozy8nK22GILhg8fTmVlJe3ateOdd96hX79+nH76\n6VxyySXstttu7LzzzixevJhJkyZx3HHHpRwfhGRRWVnJ6aefvvrYlClTGDRoEIsWLWKjpsanRzp+\n9TXjF88CIJP/LYk7WWwFzKn3fi4hgYg0ySzs4rfDDnDWWXDUBsCSxufNarcxY3vfQG0tqx81NYlf\nJ3uf7Nz0MjPrAAALvklEQVT1v9wGamc1iuFV68n2W/1nzfmrakPiqKnBalaFZLNqFfZ9vfc1q6C2\nhvasooPVsE67VbS3GjqWrKJDyZrnDtbguSQ6z8JzhwbP4bGK71Ymvu9d873z4O0LKaGGdtSGZ6+h\nhFraUUMJ0WuvoR01WPS6JDq3pO61h+tWfx69rv95SM619RJ37Zpk7bUJ48vE0NbWllFSUsKECRP4\n7W9/S48ePSgpKeGEE06gX79+HH300SxbtozBgwcze/ZsNthgAwYMGLA6WaQyz2LKlCnMnj2bM844\ng671VuccOHAg22+/Pffddx9nnHFGwmu/PPF8KClhs+WLWvUzNiXuZJGyysrK1a8rKiqoqKiILRbJ\nL039G+zUKQyiyoWjNwQWNz7epUvob1mjJHp0aHxyA+5rklL951SPNfXZ8z+9GZZ926i+rzptxAa3\nX0Ntbai77nllLXxf730unvlD/i710r17dx599NGEn5144omceOKJCT+rqalZ/fruu+9OeM7ee+/N\n8iZmr37wwQfNxvXwy5sy838zmLJyRRoDgVMX62goM9sHqHT3Q6L3wwFv2Mmt0VDSnP222IMOX3Vu\ndHzl5kt5dd47bSaGVB294Tarb1OsdXyDnoz/9j+5DygBrQ2VnvrfV/0/X4OiGQ31FtDLzHoC84DB\nQOt6h6TN2WXgoISzcndJ4152McSQqm86daVi8TaNjq/stDT3wUjBiDVZuHuNmZ0FPMuaobOfxBmT\nFJ5sL9NQKDGkqpASm6Tv+8035ei6NwlakC2lSXkiknd0Gyo9WqJcRETygpKFiIgkpWQhIiJJKVmI\niEhSShYiIjmgbVVFRCSpVLdVHTNmDO3bt6dLly506dKFXr16cfvtt2c5uuTinpQnIpKSTCytHufy\n7OkMBe7Xrx8vvfQSANOmTePHP/4x++67L7vttlu2wktKyUJECkL19OlUJtiMojLHZeRiW9X6+vTp\nw0477cQnn3wSa7LQbSgRkRTlalvV+t566y1mzJjBnnvu2eqyWkPJQkQK2+TJYenhVB6t2SYPePPN\nN5k3bx7XXnstpaWldOzYkX79+gEwcuRILrzwQsrKyigpKWH48OFMmzaNOXPmJCm1sddff52NN96Y\nLl26sM8++3DSSSfRq1evVsXeWkoWIlLY+vcPa5un8ujfv1VVpbKt6sYbb8zGG29M165dV2+rmq59\n992XhQsXsmTJEr766is+/PBDLrroolbF3lpKFiIiKaq/rWpDPXr0YOTIkSxcuJCFCxeyaNEivvvu\nO/bZZ59W1dmtWzeOPfZYJkyY0KpyWksd3CJSEErLyhJ2RJeWleWsjFxsqwprj5xasGABjz76KLvu\numva5WSSkoWIFIRMDG1tbRnZ3la1zhtvvEGXLl0AWHfddTnwwAO56aabWhV7a2mJchHJO1qiPD1a\nolxERPKCkoWIiCSlZCEiIkkpWYiISFJKFiIikpSShYiIJKV5FiKSd3r27JnWvIS2rmfPnlmvI7Z5\nFmZ2HGFl4J2Avdx9ajPnap6FiEiaimWexQfAMUDrloEsAFVVVXGH0CqFHH8hxw6KP26FHn8mxZYs\n3P1f7j4DKPq2ZqH/hSvk+As5dlD8cSv0+DNJHdwiIpJUVju4zWwSsFn9Q4ADF7l7vOvtiohIymJf\nSNDMXgTOS9bBncOQRESKRqY6uPNl6GyzP0ymflgREWmZ2PoszOxoM5sD7AM8YWZPxxWLiIg0L/bb\nUCIikv/yejSUmR1iZp+a2XQzuyDueJpiZv8xs/fM7F0zezM6tpGZPWtm/zKzZ8xsg3rnX2hmM8zs\nEzM7KIZ47zKz+Wb2fr1jacdrZn3N7P3ozydn23g1Ef+lZjbXzKZGj0PyMX4z625mL5jZR2b2gZmd\nHR0viO8/Qfy/jY4Xyve/jplNif6tfmBml0bHC+X7byr+7H//7p6XD0Iimwn0BDoA04Ad446riVg/\nAzZqcOwa4PfR6wuAEdHrnYF3Cf1F20Q/o+U43h8BfYD3WxMvMIUw+x7gKeDgGOO/FDg3wbk75VP8\nwOZAn+j1+sC/gB0L5ftvJv6C+P6jutaNntsBbwDlhfL9NxN/1r//fG5ZlAMz3H2Wu68E7geOijmm\nphiNW2lHAWOi12OAo6PXRwL3u/sqd/8PMIPws+aMu78CLGpwOK14zWxzoLO7vxWdN7beNVnVRPyQ\neKDEUeRR/O7+lbtPi15/B3wCdKdAvv8m4t8q+jjvv38Ad18evVyH8EvUKZDvH5qMH7L8/edzstgK\nmFPv/VzW/KXMNw5MMrO3zOxX0bHN3H0+hH9gwKbR8YY/1xfkx8+1aZrxbkX4M6mTD38+Z5nZNDO7\ns95thLyN38y2IbSQ3iD9vy/5FP+U6FBBfP9mVmJm7wJfAZOiX5gF8/03ET9k+fvP52RRSPZz977A\nYcCZZrY/a7J9nUIbSVBo8d4GbOfufQj/iG6IOZ5mmdn6wDhgWPQ/9IL6+5Ig/oL5/t291t13J7To\nys1sFwro+08Q/87k4PvP52TxBdCj3vvu0bG84+7zoudvgPGE20rzzWwzgKjJ93V0+hfA1vUuz5ef\nK9148+rncPdvPLr5CtzBmlt7eRe/mbUn/KL9u7s/Fh0umO8/UfyF9P3XcfclQBVwCAX0/depH38u\nvv98ThZvAb3MrKeZdQQGA4/HHFMjZrZu9L8szGw94CDCirqPA6dGp50C1P1SeBwYbGYdzWxboBfw\nZk6DDoy173GmFW/UVF9sZuVmZsDJ9a7JhbXij/6B1/n/gA+j1/kY/2jgY3e/ud6xQvr+G8VfKN+/\nmW1Sd4vGzDoBAwj9LgXx/TcR/6c5+f5z0Xvfil7/QwijLWYAw+OOp4kYtyWM1HqXkCSGR8c3Bp6L\n4n8W2LDeNRcSRiV8AhwUQ8z/BL4EVgCzgSHARunGC+wR/cwzgJtjjn8s8H70ZzGecA867+IH9gNq\n6v2dmRr9PU/770uexV8o3/8Po5inRfFeFB0vlO+/qfiz/v1rUp6IiCSVz7ehREQkTyhZiIhIUkoW\nIiKSlJKFiIgkpWQhIiJJKVmIiEhSShYiIpKUkoUUNDOridbvfzd6/n10/EUz65vhuoaZWWm990+Y\nWZcMlHu3mX1mZkPTvO4FM1ua6Z9TJJF82YNbpKWWeVjEMRfOAf4OVAO4+xEZLPt37v5IOhe4+wFm\n9kIGYxBpkloWUugSreG/9glmA8zsNTN728weiNbzOtjMHqx3Tn8zezx6fZuZvdlgJ7LfAlsCL5rZ\n89Gxz81s4+j1udH575vZsOhYTzP72MxGmdmHZjbRzNZJId67oxheN7OZUWx3RWWNTvfnF8kEJQsp\ndJ0a3IYaVP9DM+sKXAz81N33BN4BziWsA1QeLcYG8HPCBlsAf3D3cmA3oMLMdnX3vxBW5axw959G\n53lUR1/C4nN7AfsCp5nZbtE5vYC/uPuuwGLg2BR/rg3dfd8o1seBG9x9Z6C3mfVOsQyRjNFtKCl0\ny5PchtqHsDXmq9Hqmh2A19y9xswmAgPN7GHgcOD86JrBZnYa4d/H5tH1H9J4pd46PwIedfdqADN7\nBNgfmAB87u4fROe9Q9jaMhUToucPgK/c/ePo/UdRGe8nukgkW5QspNgZ8Ky7n5jgsweAswhbtL7l\n7sss7P52HrCHuy8xs7uB0gTXpmpFvdc1aZRVd11tgzJq0b9biYFuQ0mhS3bP/g1gPzP7Aazef2T7\n6LPJQF/gNNbcguoCfAcsjTbDObReWUuizxvW/TJwtJmVRnuaHBMdSyW+VKhfQmKn/6FIoSs1s6mE\nX6gOTHT3P0Svcff/mtmpwH1R57IT+jBmuHutmT1B6G84OTr/fTObRlj7fw7wSr267gAmmtkXUb9F\nXR3vmtk9hA27HBjl7u+ZWU9atj1nc1t8ak8BiYX2sxCJWXSr6wl3f7gF174InOfuUzMfmcgaug0l\nEr/FwOUtmZRH2KlxZVaiEqlHLQsREUlKLQsREUlKyUJERJJSshARkaSULEREJCklCxERSer/AQYb\nqdVnzGeaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f65d510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(elev_mec,ab['cellA'],'bs-',label='cell A')\n",
    "plt.plot(elev_mec,ab['cellB'],'rs-',label='cell B')\n",
    "plt.title('Ablation')\n",
    "plt.xlabel('Elevation [m]')\n",
    "plt.ylabel('ab [m/yr]')\n",
    "plt.legend(loc='lower right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying a linear regression on a variable that has two such distinct regimes is tricky. We better separate the two regimes before regression. The ELA (SMB==0) seems like a natural regime divider:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10fb25110>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEZCAYAAABmTgnDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOX5//H3HUCiKCgoCgoIahTFjVa+iCLRqlhx119B\nrHWhouKCdalYpUSrFdeCVluxLlj3tQqtKC6h7ihIQRCCC6uIIAiCBiTcvz+eA8aQZCYwM2dm8nld\nV65Mzpw5534GnXue3dwdERGR2hTEHYCIiGQ/JQsREUlIyUJERBJSshARkYSULEREJCElCxERSUjJ\nQrKGmR1sZh/HHUeuMrO+ZjYm7jgkP5nmWUimmdnnQD93fy3uWEQkOapZSL1nZg1ScU6q7ymSTZQs\nJGuYWQ8zm1vp78/N7DIz+5+ZLTWzx8xss0rPH2NmH0bPvWlme1d67koz+8TMlpvZR2Z2QqXnzojO\nv93MFgNDqolliJk9ZWb/NLNvgDMsGBRdd5GZPW5mW1d6zW/MbFb03DVR/IdtzPXMrHF07uKofO+Z\n2XbRc2ea2adR2T41s1MrleuNSvF0M7PxlV5/YKXnXjez66L3YbmZjTGz5pv4Tyh5TMlCsk3VdtH/\nBxwJtAf2Bc4EMLP9gfuAc4DmwD3AC2bWKHrdJ8BB7t4UuBZ42My2r3Td/4vOaQncUEMsxwFPuvvW\nwCPAxdGx7kBrYClwdxTPnsBdwKlAK6BZdM5GXQ84A2gK7BiV7zzgezPbAhgO9IzK1g2YVPX9M7Nt\ngNHAMKAF8Bfg39HxdU6N7rMd0Bi4vIb3QUTJQrLecHdf6O7fAKOA/aLj5wB/d/cPPPgnsAroCuDu\nz7j7wujxU8BMoEul685397vdfa27r6rh3u+4+6joGquAc4Gr3X2Bu/8AXAecYmYFwMnAC+7+jruv\nAf64idf7gfAhXxSV70N3XxFdpwLY28wKo/emukEBvYAyd380KuPjwHTg2ErnPODun0axPFnpvRXZ\ngJKFZLuFlR5/B2wZPW4HXGZmS6KfpcBORN/moyahdU1US4G9gG0rXWsuiVU9px3w3Lp7AtMIH+rb\nR/ddf767fw98vQnX+yfwEvC4mc0zs6Fm1sDdvwN6A+cDC8xslJntXk3srYHZVY7NJtRU1vmy0uPK\n763IBpQsJFfNBW5w9+bRzzbuvqW7P2FmbYERwIDo+DbAVMAqvT6ZYYBVz5kD/LLKPZu4+wJgASFZ\nAWBmmxNqBht1PXdf4+5/cve9CE1NxwK/AXD3se5+JLADMCMqa1VfADtXOdYWmJ9EuUU2oGQhcdks\n6sRd91PX0UH3AueZWRcAM2tiZkebWROgCbAWWGxmBWZ2FtApBTHfA/w5SkaY2XZmdlz03NPAsWbW\nNeo3KdmU65lZsZl1ipqkVhBqHGvNrKWZHRf1XfwQPbe2mmv/B9jNzPqYWQMz6w10JDTlidSZkoXE\n5d+Epo/vo98bjEiilm//7j6B0G/x16gJp4zQWUvUhn8b8C6hqWUv4M0UxDwceB542cyWAW8T9YO4\n+zTgIuAJwrf65cBXhH6UOl+PUGt4GlhGqBW9TmiaKgAuJdQQFgOHEJqkfsLdlwDHEDqtF0e/e7n7\n0nWn1Ln0Uq/FPinPzJoB/yB881sLnO3u78UalMgmimo43wC7unvVvgORnJMNNYvhwH/cvSNhaKSW\ne5CcZGHex+ZRorgNmKxEIfki1mRhZk2B7u7+AEDUqbc8zphENsHxhCaoecAuQJ94wxFJnVibocxs\nX8JIjmmEWsUHwMBo2KGIiGSJuJuhGgKdgbvcvTOho3NQvCGJiEhVDWO+/zxgrrt/EP39NHBl1ZPM\nTCM3REQ2grtb4rMSi7VmES3HMNfMiqJDvyA0SW147r/+ha9di7vn3c+QIUNij0HlU9lUvvz7SaW4\nm6EgLKb2iJlNIvRb/Lnas66+Gg49FCZOzGRsIiJCFiQLd/+fux/g7vu5+0nuvqzaEydNgr59oVcv\nOOMMWLAgw5GKiNRfsSeLpDVsCP37w4wZ0KYNrF4dd0QpU1xcHHcIaZXP5cvnsoHKJz+KfQZ3MszM\ncyFOEZFsYmZ4ijq44x4NlVrffw+bbx53FCKSYjvvvDOzZ2syfE3atWvHrFmz0nqP/KpZ9OsHCxfC\nLbdAx47pD0xEMiL6hhx3GFmrpvcnlTWL3OmzSMbdd4cRU4ccAhdcAIsWxR2RiEheyK+axTpffw3X\nXQePPAJDhsBFFyV8ydD+/SkvK9vgeGFREYNGVLe3jIhkimoWtctEzSK/+izWadEChg8PtYspU5J6\nSXlZGSXjxm1wvCTFoYmI5KL8aoaqqqgITj457ihERDZQUFDAZ599BsBZZ53FH//4x5gjql1+J4ua\nuMO8eXFHISL1mFntrUOzZs2iQYMGXHDBBRmKqHb52QyVSFkZdOsWRk9dfTU0axZ3RCJSR/37D6Ws\nrHyD40VFhYwYkdzi1am4xsZK1Afz0EMP0bx5c5544gmGDRtGo0aN0hpPIvUzWey+e+jLGDw4PB48\nmIK11e15LyLZqqysnHHjSqp5prpj6bsGwLx58xg4cCBvvPEG7s6pp57KHXfcwf3338+tt97KwoUL\n6dKlC/fccw9t27ZN6poPPfQQ119/PSUlJYwaNYqTTjqpTjGlWv1MFgCtW8N998H//geXXcaAKVMY\n3qULS6tM6issKqrhAiIisHbtWo455hgOP/xwHn74YRo0aMAHH3zACy+8wNChQxk9ejS77rorQ4cO\n5dRTT+Wtt95KeM033niD+fPn06dPH6ZOncrIkSOVLGK3774wdizbvvwyAw89FDbbLO6IRCSHjB8/\nngULFnDzzTdTUBC6gbt168bRRx/NVVddRVH0hXPQoEHccMMNzJ07lzZt2tR6zYceeoijjz6aZs2a\n0bdvX3r06MHixYvZdttt016emtTPDu6qzKBnTyUKEamzuXPn0q5du/WJYp3Zs2czcOBAmjdvTvPm\nzWnRogVmxvz582u9Xnl5OU899RR9+/YFoGvXrrRp04ZHH300bWVIhpJFIpMmwcqVcUchIlmqTZs2\nzJkzh7VV+j3btm3LPffcw5IlS1iyZAlLly5lxYoVdO3atdbrPfvssyxfvpwBAwbQqlUrWrVqxRdf\nfMHIkSPTWYyE1AyVyIMPwtNPww03wOmnQ4Hyq0g2KCoqpLqO6HA8c9fo0qULrVq1YtCgQZSUlNCg\nQQMmTJjAueeey+DBg9l3333Zc889WbZsGWPHjuWUU06p9XojR46kX79+3HDDDeuPzZs3jy5dujB1\n6lT22muvpGNLqbi3/Utya0CP1dtvu3ft6r7//u6vvRZvLCL1UOyfAQnMnTvXTzjhBG/RooVvt912\nPnDgQHd3f/jhh33vvff2Zs2aedu2bb1fv37rX1NQUOCffvqpu7ufeeaZPnjwYJ8/f743atTIp06d\nusE9evXq5VdccUW196/p/YmOp+RzOD/XhkoHd3jySRg0KMwKv/XWeOMRqUe0NlTtMrE2lJJFXZWX\nwxdfQIcOcUciUm8oWdROySKSVclCRDJOyaJ22s8ilyxbBs88E5qrRETyjJJFqixYAH/6E3TvDu+9\nF3c0IiIppWSRKnvsARMmwNlnw4knQt++oD2DRSRPxJ4szGyWmf3PzD40s/Fxx7NJGjQIyaKsDHbb\nDTp31lLoIpIXYu/gNrPPgJ+5+9JazsnNDu4lS6B587ijEMl56uCuXX3p4DayI47UU6IQkTyRDR/S\nDow1s/fN7Jy4g8mIv/0t6b3BRSQ/5dq2qtmwNtRB7r7AzLYjJI2P3f3NqieVlJSsf1xcXExxcXHm\nIkw1d/jFL+CEE+C662CHHeKOSEQyrKZtVceNG8dhhx1GkyZNAGjWrBn9+vX7yWdgTUpLSyktLU1h\nlD+KPVm4+4Lo9yIzew7oAtSaLHLegAFw6qlhccJOneB3v4NLL4UqGy+JSM2G9u9PeVnZBscLi4oY\nNGJExq6xsWrrg9lxxx2ZM2cOEJY6P/jgg+ncuTPHHXdcrdes+kX62muvTUmsEHOyMLMtgAJ3X2Fm\nTYAjgdSVLptts01YX+r888N6U9dcA7fdFndUIjmjvKyMknHjNjhekuFrQHq2VV2nXbt2dOvWjWnT\npiVMFukUd5/F9sCbZvYh8C4wyt1fjjmmzNplF3jqKbjpprgjEZGNsG5b1fbt2zN79uz126Gu21b1\nX//6F4sWLaJ79+6ceuqpdb7+zJkzeeuttzjwwAPTEH3yYq1ZuPvnwH5xxpA1GsbeIigiGyEd26rO\nnz+f5s2bU1FRwYoVKzjxxBM56KCD0l6W2sRds5DaTJwY+jKW1jgFRURiluptVSH0WSxZsoRly5bx\nzTffUFhYyG9+85t0FSEpShbZbKed4LvvYPfdYfhwWL067ohEpIpUb6ta1VZbbUXfvn0ZPXp0KsOu\nM7V9ZLOWLeHvf4cLL4QrroC77oKbb4bjj4caht2J1BeFRUXVdkQXRs0+mbpGqrdVhZ+OlFqxYgWP\nPfYYnTp1SjqmdFCyyAWdOsGLL8JLL8H110OPHmE0lUg9loqhram4RkFBAaNGjeKiiy6ibdu2FBQU\n0LdvX4YNG8bKlSvp06cPc+bMoVmzZhxxxBHrk0VN8ywAFixYQNOmTQFo3LgxXbt25eGHH97kWDdF\n7GtDJSNn14YSkZTQ2lC1qy9rQ4mISJZTssgH7mHpkPvug4qKuKMRkTykZJEPzODqq2HkyLCHxtix\ncUckInlGfRb5xB2eew5+//sw3PaWW2DPPeOOSmSTqc+idpnos1CyyEerV8Pdd4f9NGKeyCOSCkoW\ntVOyiChZiNRvSha1y0Sy0DyL+qaiAgoKNKlPckq7du1qnZdQ37Vr1y7t91DNop5Yt25/p6++ouu8\neby0yy7MbdYsI+v2i0g8VLOQOqu6bn+/SZMAGP7993GFJCI5RENn67lzJk4M606Vl8cdiohkMSWL\neu7un/8cmjaFzTaLOxQRyWJqhqrnVjRuDIMHxx2GiGQ51SykZtp0SUQiqlnUE3Vet7+iArp1gwMP\nDMuit26dzvBEJMtp6KzUbNkyuPFGuPdeuPhiuPxyaNIk7qhEJElaolwyo1kzGDoUJkyA6dPDelNj\nxsQdlYjEQDULSd5774WaRczbO4pIcrQ2lIiIJJR3zVBmVmBmE83shZrO+fLLTEYkdbJoESxeHHcU\nIpJGWZEsgIHAtNpO6NgRzjgDJk7MUESSvLFjwz/QrbfCqlVxRyMiaRB7sjCznYCjgX/Udt6nn8Je\ne8Hxx8Mhh0ClZY4kbn37whtvwH//GzZbeuqpsBGTiOSN2PsszOwp4AagGXCZux9XzTnr+yx++CFs\nBteyJRQXZzRUScZrr8Fll4UlRF59FRpqKo9IXPJm1Vkz6wUsdPdJZlYM1FiokpKS9Y+Li4spVqbI\nTocdBh98AOPHK1GIZFhpaSmlpaVpuXasNQsz+zPwa2ANsDmwFfCsu/+mynlJjYZatgz69YPzzw+f\nWdorRUTqs7wZDeXuf3D3tu7eAegDvFY1UdRF48bQsycMHAj77AP33QfariHLvPQSrFkTdxQiUkex\n91msY2Y9SKLPIhnuobl82LDQGnLnndC7d+2v6d9/KGVlG+7pUFRUyIgRg5K+t9Ri1So4+ugwDvqW\nW+CXv1T1TySNNCmvDsrKwu+a1stbp7i4hHHjSjY43qNHCaWlGx6XjeQOo0eHDZfatg3DbffZJ+6o\nRPJS3jRDZUJRUc2JIgfyZP4xg2OPhSlTwjjoI46Ap5+OOyoRSSDvk0VNFi6EXXeFm2/Wtg2xaNQI\nLrgAZsyAI4+MOxoRSaDeJovtt4cnn4SPPoIOHWDmzF5xh1Q/bb11mJMhIlmt3iYLgJ/9DB56CKZN\ng4YNv4s7HKnsnXc0TV8ki2jWFNCqFfTs+Q7Tp79BRUUjGjX6cbxtUVFhjJHVY8uWwXnnwf77h7bC\n3XaLOyKRei3vR0OlyldfhSVGJIPKy2H48DDM9vTTYfBgaN487qhEcoZGQ2XYl1+GRVX79IF33407\nmnqksBCuvDK0E65aBYcfriFsIjFRzSJJy5fDAw/AHXfAdtuFWeKnnBIG9UiGfPcdbLFF3FGI5AxN\nyotRRUWYUzZsWJiMfMUVcUckIlI9JYsssXYtFKghL16rV8NVV8Ell0CbNnFHI5JV1GeRJapLFOvW\npVq7NvPx1EsVFdCkCey3H1xzDXz7bdwRieQlJYsU++qr0DTVsSPcfTesWBF3RHlu883huutg0iSY\nMwd23x3uvTckERFJGSWLFNt+e5gwIXxevfoq7LxzSB5z5sQdWZ5r0ybMsBw1Ch55JCyFLiIpoz6L\nNJs1C/76V9h33zBVQDJg3X8rWv5c6jl1cIuISELq4M4Tq1bBTTeFfg7JgGHD4MYbtX2iyEbQ2lAx\nWrECPvkk9MmedNKP28FKmhxzTJgR3rEj/PnPDH39dcpnztzgtMKiIgaNGBFDgCLZS81QWWDRIhgx\nAu66C/bYIwzuOfjguKPKY2+8AZdeyrwZM9ipmqG2JT16UFJamvm4RFJMzVB5Zrvt4OqrQ2d4v34a\n9Zl23bvDe+/x3o47xh2JSM5Qssgim20Gp50GPXrEHUk9UFDAlO23jzsKkZyhZJEjli6F3r3DfkB5\n3CKXPfQmi/yEkkWOKCyEQw+Fc8+Fzp3hwQfDdg+SeluuXh02XXr+eSUNkUiNHdxmdkcSr1/u7tek\nNqRqY8nrDu66WLsWXn45jAKdNCl0ip98ctxR5aah/ftTXla2wfHCoiIGnXQSXH556FC67baQoUVy\nTEYm5ZnZbOCPCV4/yN07bvTNzRoD/wU2Iwzjfdrdr63mPCWLaoS9w6GoKO5I8tSaNXD//TBkCPTs\nGeZotGoVd1QiSUtlsqhtnsVf3H1kgkC22ZSbu/sqMzvU3b8zswbAW2b2oruP35Tr1hd77lnzc+5a\n7WKTNWwI/fuHLRKHDg2bL4nUU7XOs4g+wC9297+kPRCzLQi1jPPd/f0qz6lmUQdffBFGVA0YAGef\nDc2axR2RiMQhY/Ms3L0C6JuKG9XEzArM7EPgS2Bs1UQhdde6NTz8MLz/PrRvH2aGf/JJ3FHlKY0y\nkHoimeU+3jSzvwJPACvXHXT3iakIwN3XAvubWVPgX2a2p7tPq3peSUnJ+sfFxcUUFxen4vZ56//+\nDx59FObNC/tqdOsG118fWlUkhQYOhPnz4ZZbwjIiIjEqLS2lNE2rDyRc7sPMXq/msLv7YSkPxmww\nsNLdb69yXM1Qm+i778KX4ObN444kz6xaFdagHzo0TIQZMiSMoBLJAple7uNwdz+0yk9KEoWZbWtm\nzaLHmwNHANNTcW35qS22qDlRLF6c2VjySuPGcNll8PHHYZ/djh3hzjvjjkok5ZKpWXwGPAPc7+4f\np/TmZnsDIwlJqwB4wt1vqOY81Sw2Uf/+Qykr27B9fccdt+U//7mQXr3gkkvg5z+PIbh8MmNGmADT\nu3fckYhkbOjsOvsCfYD7zKwAuB943N2Xb+rN3X0KoNlOGVBWVs64cSUbHO/Ro4TPPoP77guT+9q0\nCUnjhBPCyFGpo913Dz8ieSZhM5S7f+vu97p7N+BKYAiwwMxGmtmuaY9Q0m6bbcJk5U8/DYli2DC4\nI5n5+5I89zCmWSRHJUwWZtbAzI4zs+eAYcBtQAdgFPCfNMcnGdSwIZxyCrz5ZhjkIyn06aew995h\n86Vly+KORqTOkungngkcD9zi7vu7++3uvtDdnwbGpDc8iUuDBhseW7sWXn9da+ttlF13hSlTwk5X\nRUVhPPOaNXFHJZK0GpOFmZ1qZi2Afdy9n7u/XfUcd784rdFJVvnyy9BM1alT2NlPq1/UUevWYa2p\nl16CZ54JNQ3NlpQcUVsXZlvgKaCRmb0KvAiM17Ck3FRUVAiU1HA8Oa1bh4E+paWhX+Pqq+G3v4UL\nLoCddkpZqPlvv/3glVdgzBi9cZIzkhk6uxVwOHAU0AX4mND89JK7L0x7hGjobLb65JMwpeCAA+DX\nv447GhGpKiNLlNdy8z2BXwJHunvPVASRxD2VLKT+mDIFOnSAJk3ijkRyXKZncGNm+0Qjok4C9gA+\nz1SikNxUXh72DPr667gjyUEPPgh77AEjR4ZRBSJZIJmhs/cTJuKdDBwb/RyT5rgkx337bfiCvOuu\ncN55YaMmSdJtt8ETT8Df/x6m1L9e3fJsIpmVTJ/FNHevZZud9FMzVO5auDB85v3tb7DvvnDttdC1\na9xR5Qh3ePJJGDQITjwRbr898WtEKslon4WZ3QfcVt2y4ZmiZJH7Vq0KX5bbtQsbM0kdlJfD3Lmw\n225xRyI5JtPJogfwAmFzolWAEZYo3ycVASRDyUJEpO4y3cF9H3A6Yejsuv6KY1NxcxGAJUugb194\n6y3NDq+Tb7+FZ5/VmyYZkUyyWOTuL7j75+4+e91P2iOTeqNxYzjwQDjjDOjSBR55BFavjjuqHLBg\nQegE6t4dxo+POxrJc8k0Q90NbE1YOHDVuuPu/mx6Q/tJDGqGqgcqKuA//wmzw6dPD8snHX983FFl\nuYqKMMT2mmvg0EPhxhuhbdu4o5Iskek+iweqOezufnYqAkiGkkX9M3kyFBaGNfckCStWwM03w113\nwcSJYSSB1HuxzuCOg5KFVOYOlpL//PPQ119DixZxRyFZIiMd3GbWP4lAEp4jkkrz5v24zfW338Yd\nTRZSopA0qbFmEe29fXltrwWuc/e90hFYlVhUsxAg1CreeSf0a7z6augUv+giaN8+7siy3D/+EWZD\nduoUdySSQZkaOjuOH5f3qO7nGGBsKoIQSZYZdOsWJjZPnBh29zvggLCHuNRi9Wo47DA499wwrV6k\njtRnITlvxYrwWdi8edyRZLklS+D668Poqcsug9/9DjbfPO6oJI3UwS2SpCVLlEQ28MknYb2pHXeE\n4cPjjkbSSMlCJAlz58I++8AJJ8DAgWGDusr69x9KWVn5Bq8rKipkxIhBGYoyRj/8AI0axR2FpFEq\nk0Vt26qmnZntBDwEbA+sBe519zvijEnyR5s24Uv0vffCMceEdfgGDoRjj4UGDaCsrJxx40qqeWV1\nx/KQEoXUQTL7WbQwszvNbKKZTTCz4WaWqvF5a4BLoxFVBwIXmNkeKbq2CC1ahBaXzz8Pfbs33RTm\nrUkNJk8O/RlLl8YdiWSZZNaGehz4irD50SnAIuCJVNzc3b9090nR4xWE/b13TMW1RSpr1Aj69AnD\nbi+4IO5ostgOO8DKlWGnvjvuCE1VIiSXLFq5+5+ihQQ/d/frCc1GKWVmOwP7Ae+l+toilTVoEHcE\nWaxly7Bb1auvhoW69toLnn9eK9tKUn0WL5tZH+DJ6O9TgJdSGYSZbQk8DQyMahgbKCkpWf+4uLiY\n4uLiVIYgst6CBftTXh7Wpqq3OnWCMWPgpZfCyrYHH6zZ4TmgtLSU0tLStFy7thnc3wJOmKndhNAB\nDaE2ssLdm6YkALOGwGjgRXevdhyfRkNJOlQdDeUOS5fuwjffdGPVql3o3x8GDIBWrWIMUmQT5NXQ\nWTN7CFjs7pfWco6ShWTU9Olh/alu3eC00+KORmTjZDxZmNk2wG7A+oq5u/93k29udhDwX2AKoRbj\nwB/cfUyV85QsRLKJO/TuDT17wplnqiMoS2V6P4vfAgOBnYBJQFfgHXc/LBUBJEPJQrJJeXnoAz7z\nTNh667ijidH778Oll8Ly5XDbbXD44XFHJFVkeg/ugcABwGx3PxTYH/gmFTcXyUXLl8OECdChA1x4\nIZSVxR1RTA44AP77XxgyBM47D3r1go8/jjsqSZNkkkW5u5cDmFljd58O7J7esESyV8uW8M9/wkcf\nwTbbhIFCvXqFL9r1jhmcdBJMmxZqFm+/HXdEkibJNEM9B5wFXAIcBiwFGrn70ekPb30MaoaSrPX9\n9/Doo2E5kUMOiTsakR/FNhrKzHoAzYAx7r46FQEkeV8lC5FctXZtqIFoL9yMy3SfxXruPs7dX8hk\nohDJZV9/DaefDuPHxx1JjJ57LoxBfueduCORTVCnZCEiddO4MXTuHEaZdusGTzxRD5dbOvHEMLvx\nV78Kb8Tnn8cdkWwEJQuRNNpyy7Ah3SefwBVXwN13h1FUo0fHHVkGFRSE6tWMGWEZkZ//HH7/+9DZ\nIzkj9hncyVCfheSTDz+EJk2gqCjuSGKyYAHccw8MHqzJfGmWV8t9JEPJQuoLd/UDS+rE1sEtIukz\nezbsvXeYHb5yZdzRxGTZsrgjkBqoZiGSJdzDhOhhw+CNN6Bfv7BRU9u2cUeWnKH9+1NezXT2wqIi\nBo0YkfgCa9eGbNm1K/zpT9C6dRqirF/yZg9uEfmRGfToEX4++wz++lfYf3+4/XY444y4o0usvKyM\nknHjNjhekuwFCgrCDPAbbwxJY+DAsMVrkyapDFM2kpqhRLJQhw4hSXz+ORxzTNzRZFCzZjB0aFh8\n6+OPYffd4cUX445KUM1CJKs1rWWLsaVLw9pUeWnnneGxx+Ddd+v5loXZQzULkRw0a1aofZxzTljQ\nMG917Qr77Rd3FIKShUhO2nnnsDR6u3Zw5JFwxBHw73+HPuJ64euvw49kjJqhRHLUdtvBNdeEydBP\nPhm2lZg9O6ysEYfCoqJqO7ML0zH78JVXwmYigwaF340bp/4e8hMaOiuSJ9yhogIa1pevgNOnh0w5\ndSrcdBOcfLJmNFahGdwikrSKCnjvPTjwwDz9LH3ttTDEtmnTUONo1CjuiLKGkoWIJG3OnNCnseWW\nYepC79552GpTURGWQD/44LgjySpKFiJSJ2vXwpgxYXb4lClw/vlh2+yWLeOOTNJJa0OJSJ0UFMDR\nR8PLL8PYsTBvHrz6atxRZcgrr8CaNXFHkfNUsxCR/LV6Nfzyl/Dll3DrrXDUUXnacVO9vKpZmNl9\nZrbQzCbHHYtIfff992FNquXL444kRTbbLNQshg6FSy4JyWLKlLijykmxJwvgAaBn3EGISFgh/M03\noX378Nn66adxR5QCZnDssWGq+7HHwuGHh/1tpU6yohnKzNoBo9x9nxqeVzOUSAbNnRu2gP3HP8Le\n4UOGhL0WiZVDAAAOMklEQVTE88I334QE0qxZ3JGkXd6NhlKyEMlOK1fCww/DnntC9+5xRyN1VS/3\nsygpKVn/uLi4mOLi4thiEakvmjSBc8+NO4oMGT8eysvhkEPijmSjlZaWUlpampZrq2YhIhtl0aKw\n2sZFF4Umqv79h1JWVr7BeUVFhYwYMSiGCOvopZegf3/42c/C8iG77RZ3RJssH2sWFv2ISI5o3Bg6\ndoQTTgir4H79dXumTetdzZklGY5sI/XsGdabGj48rI1y+ukweDA0bx53ZFkh9tFQZvYo8DZQZGZz\nzOysuGMSkcSaNg01i08/DQu/zpt3YNwhbbrNNw8r2U6bFpqkiovr0brvtYs9Wbh7X3dv7e6N3b2t\nuz8Qd0wikrxGjeBXv4L9978/7lBSp2VL+NvfwnpTBbF/TGYFvQsikhE52e3YpEncEWQNJQsRSbvP\nPoN994V77w2zxHPaDz/AlVeGySj1SLZ0cItIjisqKqS6zuyiokLat4fbbw+r3l59ddg7fMAA2HHH\njIe56SoqwjIi++0XCvH738NWW8UdVdplxdDZRDR0ViR/lJXBnXeGyX5//SucdlrcEW2kuXND5nvl\nFbj2Wjj7bGjQIO6ofiLvZnAnomQhkn+++SZ8SW/RIu5INtEHH4Sd+i69FI4/Pu5ofkLJQkTylntY\n9Tanlm5a9/mUZcuf5+OkPBERIHSGH3BA2P714ovDxL+sl4YkMbR/f8rLyjY4XlhUxKARI1J+v0SU\nLEQkq+yyS5gT97e/hTlxnTuHvcOPPDIHpzzcdRd8+21Y772wsE4vLS8ro2TcuA2Ol6QotLrKtbde\nROqBHXYIfcazZ4cJf1deGYbd5pyePeH992GPPeCxx3J0skmgZCEiWauwEM46CyZNCoONcs6uu8Iz\nz8BDD8Ftt0HXrvDWW3FHtVGULEQk65mFZUWqqqiA997LgS/shxwSlkC/6CKIob8hFZQsRCRnzZkD\nffuGL+yPPRYmV2etggL49a9h5Mi4I9ko6uAWkZzVvn2Y5Dd6dFhZ/IorwqTq/v1h223jjq6O3H8y\nqqqwqKjazuzCoqKMhVSZ5lmISN743//gjjvgiCOgT5+4o6mDRYtCZ3hJCRx7bMqG4mpSnohIvhkz\nJswE33770Bm+//6bfMlUJgv1WYhIvfDdd3D33bBiRdyR1OCoo0LVqHdvOPpoOPNMWLAg7qjWU7IQ\nkXph6VJ47TVo1w4uvxxmzYo7omo0bAjnngszZkDr1rBsWdwRradmKBGpV2bNCqvdPvAAHHoo/PGP\nsM8+cUeVHuqzEBHZRN9+G0ax7rcfHHxw3NHUwapV0LhxUqcqWYiI1FcDBoS9NG65JSwjUgt1cIuI\npNFXX4Xd/CZPjjuSavzlL9CjR6gOXXghLF6ckdsqWYiIVNG4MXToEAYlHXYYvPBCWFokKzRuHHro\np08P8zH22CNMLkkzNUOJiNTghx/g6afDl/mvv4a//z1M+Msq06eH3fp+/esNnsqrPgszOwoYRqjl\n3OfuN1VzjpKFiMTGHd59NywhsttucUeTvLxJFmZWAJQBvwC+AN4H+rj79CrnKVmIiNRRPnVwdwFm\nuvtsd/8BeBzIrh3PRURq8cknYTe/Bx+E8vK4o0mfuJPFjsDcSn/Pi46JiOSEDh3gz3+Gxx+HnXcO\nawF++WXcUaVezixRXlJSsv5xcXExxcXFscUiIrJOQUFY1umoo+Djj8PApI4d4Z57wpawmVRaWkpp\naWlarh13n0VXoMTdj4r+HgR41U5u9VmISC5ZsiR0irdoEW8c+dTB3QCYQejgXgCMB05194+rnKdk\nISI5zz0sM9K0aWbulzcd3O5eAVwIvAxMBR6vmihERPLFzJmhX+Oii8LjXBJ3BzfuPsbdd3f33dx9\naNzxiIikS1ERTJkSahbduoVN8V59NdQ4sl3sk/KSoWYoEck3330HjzwCw4bBpZdCv36pv0fe9Fkk\nS8lCRPKVO6xZA40apf7aedNnISJS35lVnyjWrIEJEzIfT02ULEREstDs2XDyyXDQQfDUUyF5xEnN\nUCIiWWrNmrA8+rBhYTvYCy+E3/4WmjdP7vXqsxARqWcmTIDhw6FXL+jdO7nXKFmIiEhC6uAWEZH1\nVq6EESPCcNx0UbIQEclxS5bAv/8N7drBVVfBvHmpv4eShYhIjmvTBp5/Puzm9/33sM8+MGBAau+h\nPgsRkTyzfDlMngzdu6uDW0REElAHt4iIZJSShYiIJKRkISIiCSlZiIhIQkoWIiKSkJKFiIgkpGQh\nIiIJKVmIiEhCShYiIpKQkoWIiCQUW7Iws1PM7CMzqzCzznHFISIiicVZs5gCnAiMizGGrFBaWhp3\nCGmVz+XL57KByic/ii1ZuPsMd58JpGSRq1yW7//B5nP58rlsoPLJj9RnISIiCTVM58XNbCywfeVD\ngANXu/uodN5bRERSJ/b9LMzsdeAyd59YyznazEJEZCOkaj+LtNYs6qDWwqSqsCIisnHiHDp7gpnN\nBboCo83sxbhiERGR2sXeDCUiItkvq0dDmdlRZjbdzMrM7Mq449kYZraTmb1mZlPNbIqZXRwd38bM\nXjazGWb2kpk1q/Saq8xsppl9bGZHxhd9csyswMwmmtkL0d95UzYAM2tmZk9FMU81s//LlzKa2e+i\nybGTzewRM9ssl8tmZveZ2UIzm1zpWJ3LY2ado/ekzMyGZbocNamhfDdH8U8ys2fMrGml51JXPnfP\nyh9CIvsEaAc0AiYBe8Qd10aUYwdgv+jxlsAMYA/gJuD30fErgaHR4z2BDwn9STtH74HFXY4EZfwd\n8DDwQvR33pQtivtB4KzocUOgWT6UEWgNfAZsFv39BHBGLpcNOBjYD5hc6VidywO8BxwQPf4P0DPu\nstVSvsOBgujxUODGdJQvm2sWXYCZ7j7b3X8AHgeOjzmmOnP3L919UvR4BfAxsBOhLCOj00YCJ0SP\njwMed/c17j4LmEl4L7KSme0EHA38o9LhvCgbQPQtrbu7PwAQxb6M/CljA6CJmTUENgfmk8Nlc/c3\ngaVVDtepPGa2A7CVu78fnfdQpdfEqrryufsr7r42+vNdwucLpLh82ZwsdgTmVvp7XnQsZ5nZzoRv\nBe8C27v7QggJBWgZnVa13PPJ7nL/BbiCMH9mnXwpG0B7YLGZPRA1tY0wsy3IgzK6+xfAbcAcQpzL\n3P0V8qBsVbSsY3l2JHzerJNLnz1nE2oKkOLyZXOyyCtmtiXwNDAwqmFUHVmQcyMNzKwXsDCqOdU2\nvDnnylZJQ6AzcJe7dwZWAoPIj3+/rQnfutsRmqSamNlp5EHZEsi38gBgZlcDP7j7Y+m4fjYni/lA\n20p/7xQdyzlRFf9p4J/u/nx0eKGZbR89vwPwVXR8PtCm0suzudwHAceZ2WfAY8BhZvZP4Ms8KNs6\n84C57v5B9PczhOSRD/9+hwOfufsSd68AngO6kR9lq6yu5cm5cprZmYTm4L6VDqe0fNmcLN4HdjWz\ndma2GdAHeCHmmDbW/cA0dx9e6dgLwJnR4zOA5ysd7xONSmkP7AqMz1SgdeHuf3D3tu7egfDv85q7\nnw6MIsfLtk7UfDHXzIqiQ78AppIH/36E5qeuZlZoZkYo2zRyv2zGT2u6dSpP1FS1zMy6RO/Lbyq9\nJhv8pHxmdhShKfg4d19V6bzUli/u3v0EPf9HEUYPzQQGxR3PRpbhIKCCMJrrQ2BiVK7mwCtR+V4G\ntq70mqsIIxc+Bo6MuwxJlrMHP46Gyrey7Uv48jIJeJYwGiovyggMieKcTOj8bZTLZQMeBb4AVhGS\n4VnANnUtD/AzwjYKM4HhcZcrQflmArOjz5aJwN3pKJ8m5YmISELZ3AwlIiJZQslCREQSUrIQEZGE\nlCxERCQhJQsREUlIyUJERBJSshARkYSULCSnmVlFtMDfh9Hv30fHXzezzim+10AzK6z09+jKewds\nwnUfMLPPzKx/HV/3mpl9m+pyilQnW/bgFtlYKz0s8JcJlwD/BMoB3P2YFF77cnd/ti4vcPfDzOy1\nFMYgUiPVLCTX1bbabTjB7Agze9vMPjCzJ8xsCzPraWZPVjqnh/2409/dZjbews6GQ6JjFxFWZn3d\nzF6Njn1uZs2jx5dG5082s4HRsXZmNi1a1vwjMxtjZo2TiPeBKIZ3zOyTKLb7omvdX9fyi6SCkoXk\nus2rNEP9v8pPmlkL4BrgF+7+c2ACcClhraAuZrZ5dGpvwgZbAH9w9y6ENaGKzayTu99JWJmz2N1/\nEZ3n0T06ExaoOwA4EDjHzPaNztkVuNPdOwHLgJOTLNfW7n5gFOsLwG3uviewj5ntk+Q1RFJGzVCS\n675L0AzVlbC95FvRCpuNgLfdvcLMxgDHmtkzQC/Cyp0QVuo8h/D/xw7R6z9iw9VM1zkYeM7dywHM\n7FmgO2H13c/dfUp03gTC9pbJGBX9ngJ86e7Tor+nRteYXN2LRNJFyULynQEvu/tp1Tz3BHAhYZvK\n9919ZbSb4WXAz9x9uZk9ABRW89pkVV4yuqIO11r3urVVrrEW/X8rMVAzlOS6RG327wIHmdkuAFF/\nxW7Rc+MIGxmdw49NUE2BFcC30YY5v6x0reXR81Xv/QZwQrQvRBPgxOhYMvElQ/0SEjt9Q5FcV2hm\nEwkfqA6Mcfc/RI9x98XRLmKPRZ3LTujDmOnua81sNKG/4TfR+ZPNbBJh/f+5wJuV7nUvMMbM5kf9\nFuvu8aGZPUjY88KBEe7+PzNrx8Zt4VnbtqbaU0Biof0sRGIWNXWNdvdnNuK1rwOXufvE1Ecm8iM1\nQ4nEbxlw3cZMygPaAz+kJSqRSlSzEBGRhFSzEBGRhJQsREQkISULERFJSMlCREQSUrIQEZGE/j/D\njdmjv9I4egAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fbe53d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Only select points that experience net ablation\n",
    "from scipy.stats import linregress\n",
    "gradient = {}\n",
    "intercept = {}\n",
    "for cell,color in [('cellA','b'),('cellB','r')]:\n",
    "    x_pos_smb = smb[cell] < 0.\n",
    "    gradient[cell], intercept[cell], r_value, p_value, std_err = linregress(elev_mec[x_pos_smb],ab[cell][x_pos_smb])\n",
    "    plt.plot(elev_mec[x_pos_smb],ab[cell][x_pos_smb],color+'s',label=cell)\n",
    "    xx = np.arange(0,1201)\n",
    "    plt.plot(xx,intercept[cell]+gradient[cell]*xx,color+'--')\n",
    "\n",
    "plt.xlabel('Elevation [m]')\n",
    "plt.title('Linear regression')\n",
    "plt.ylabel('ab [m/yr]')\n",
    "plt.legend(loc='upper right')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'cellA': -0.0038208347471587302, 'cellB': -0.0049046200388104672},\n",
       " {'cellA': 4.1118660565964289, 'cellB': 5.9436234201181755})"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient, intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accumulation strangely enough also shows a elevation gradient, although each elevation class receives the same amount of precipitation. There is a simple explanation for this, namely the rain repartitioning scheme in CLM was enabled in the experiment that these data derive from. Lower elevation classes generally have higher temperatures so less rain is converted to snow. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10fe1c9d0>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEZCAYAAABSN8jfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lOXZ9/HvGSxEVFRwX8AFUVFRsaaIVmLdcEOtYtE+\norYKdUFr+6jUpYxWW3CrWvUt1KXQWtHHFVwQqgZXEEUEFRoKyia4gYhAEJLz/eO6EyZhEmYyWyb5\nfY5jjsxc9zLn3IE5c1+ruTsiIiLpKMp3ACIiUviUTEREJG1KJiIikjYlExERSZuSiYiIpE3JRERE\n0qZkIpInZvawmd2UxvErzGy3zEUk0nhKJtIsmFmZmS01sx/kO5ZsMLNXzewX8WXuvoW7f5qnkERq\nUTKRgmdmnYAjgCqgT57DEWmRlEykOegPvA38HTi/utDMis3sDjP71MyWmdlrZtYm2naEmb0Zlc8z\ns/5Rea07ADM7z8xej3tdZWYXm1m5mS03s5vMbI/oXN+Y2Wgz2yTRsXHH71H3A5jZVmY21sy+MLOv\no+c7RdtuBn4M3Gtm35rZPXXPZWbtzGxUdPwnZnZd3c9gZrdFd29zzKx3mtdcpBYlE2kO+gP/BP4F\nHG9m20bldwAHAz2A9sDVQJWZdQReAO4GtgEOAqY1cP66cw4dF3feq4HhwDnArsABwNkNHFvf/EVF\nwEPROToCq4D7ANz9euB14DJ3b+fulyc4173AFsBuQCnQ38wuiNteAswEOgC3AQ/W92FFGkPJRAqa\nmR1B+PJ93N2nAv8FzjEzAy4ALnf3JR5Mcve1hC/+Ce7+uLtXuvsyd5+ewtsOc/eV7j4T+BAY7+7z\n3H0F8CIh0dQbcqJCd1/q7k+7+xp3Xwn8CThyYx8fwMyKgJ8Bg919lbvPIyTSc+P2nefuD3mYjG8k\nsIOZbZfEZxVJipKJFLr+hC/zZdHrR4HzCHccxcDcBMfsCsxJ4z2/iHu+Gvi8zuvNUz2hmW1qZsOj\nKrlvgInAVlFS3JhtgE2A+XFl84Cd414vqX7i7qsJiSjlOEXqs0m+AxBpLDMrBs4CisxscVTcBtgS\n2JHwxb4nMKPOoQsI1T6JrATaxr3eIY0Qa53LzBo61/8CewGHuvuXZnYgMJXwpe/UXz0G8BWwFugE\nzIrKOgGLGh+6SGp0ZyKF7HRgHbAvcGD02IfQvtCf0AbxZzPb0cyKzKxH1HX4EeBoMzvTzFqZWfvo\nyxtC28lPozuFzsAv04jvA2A/M+sWNfwPof6ksDkh+X1rZu2BWJ3tnwMbNNwDuHsV8Dhwi5ltHvVu\nuxL4Rxqxi6REyUQKWX/gIXdf5O5fVD8IDdfnAIMJdyVTgK+BoUCRuy8ATiTcDSwF3ge6Ref8M+Gv\n/CXAw4SG/XjJNqjj7rOBm4CXgXJCkqvPXYS7mK+AtwgdBOLdDfSNenrdleC9Lyc02s8FXgP+6e4P\nN/B+WshIMsqyvThW1AXxLkLietDdh9XZ3gt4lvV120+5+83JHCsiIk1DVpNJ1MukHDga+IzwF2I/\nd58Vt08v4Lfu3ifVY0VEpGnIdjVXCTA76ja5FhgNnJpgv0Q9VpI9VkRE8izbyWRnQs+Zagup3V2x\n2mFmNs3MnjezrikeKyIiedYUuga/B3R091VmdgLwDNAlzzGJiEgKsp1MFhFGJ1fbhTp93939u7jn\nL5rZ/VHXyI0eW83M1DNFRCRF7p7MoNikZLuaawrQ2cw6mVlroB8wJn4HM9s+7nkJoVPA0mSOjefu\nBfkYMmRI3mNQ/PmPQ/EX5qOQ48+0rN6ZuHulmV0GjGd9996ZZjYwbPYRwJlmdjGhb/9qwhxD9R6b\nzXhFRKRxst5m4u7jgL3rlA2Pe34f0eyoyRwrIiJNj0bA51lpaWm+Q0iL4s8vxZ9fhR5/JmV9BHwu\nmJk3h88hIpIrZoYXUAO8iIi0AEomIiKSNiUTERFJm5KJiIikTclERETSpmQiIiJpUzIREZG0KZmI\niEjalExERCRtSiYiIpI2JRMREUmbkomIiKRNyURERNKmZCIiImlTMhERkbQpmYiISNqUTEREJG1K\nJiIikrZN8h2AiDQtQwcMoKK8fIPy4i5dGDxiRB4iSk2hx59N9V2bTFAyEZFaKsrLiU2cuEF5LPeh\nNEqhx59N8dfmxgyfW8lEJAea3F/LVVXw3XewYsX6x7ffhp9LliQ+Zt48uOWW3MbZGPPm1V9eCPFn\nU33XJgOynkzMrDdwF6F95kF3H1bPfocCbwE/c/enorJPgeVAFbDW3UuyHa9INmTkr+V16xJ/+dd9\n3tC26uerVsGmm0K7drDFFuFR/XzZssTvX1kZjmvqKivrLy+E+LOpvmuTAVlNJmZWBNwLHA18Bkwx\ns2fdfVaC/YYCL9U5RRVQ6u71/OsWKXALFsBNNyWXGNasWf/FH//lX/f1llvCrrsm3q/6+WabQatW\niWMqLYUvvtiwfI89CuMv+zffDNe1rkKJP5vquzYZkO07kxJgtrvPAzCz0cCpwKw6+w0CngAOrVNu\nqMeZFLLvv4eyMpg9O/H2NWtg7VrYbjvo3LnhBNC2LZjlNHyRZGU7mewMxKfBhYQEU8PMdgJOc/ej\nzKxuNZYDE8ysEhjh7n/LarQimbB0Kbz4Ijz7LIwfD127Qps2ifft3Bn+8IfcxrcRxV26JKx+K+7S\nJdehNEqhx59Nta5NgmrXdDSFBvi7gGviXsf/6XW4uy82s20JSWWmu7+R6CSxWKzmeWlpKaWlpVkI\nVaQec+eG5DFmDLz3Hhx1FPTpA3/5C2y/fag6+uSTfEeZlELvPlvo8WdLWVkZFTvtBDvtFAoynEzM\n3TN6wlonN+sBxNy9d/R6MODxjfBmNrf6KbANsBIY4O5j6pxrCLDC3e9M8D6ezc8hsoGqKnjnnZA8\nxoyBr76CU04JCeToo0OVVJwm15tLWjwzw90zVm+a7WTSCvgPoQF+MfAOcLa7z6xn/4eBse7+lJm1\nBYrc/Tsz2wwYD9zo7uMTHKdkItm3ejX8+98heYwdCx06wKmnhgRSUgJFat6TwpHpZJLVai53rzSz\nywiJoLpr8EwzGxg2e90/yeIzwvbA02bmUZyPJEokIln1+efw/PMhgbz6KnTvHpLHNdeE9g4RAbJ8\nZ5IrujORjHGHWbNC8nj2Wfj4YzjuuJBATjwR2rfPd4QiGVFQ1Vy5omQi9UmqrWLdOnjrrfUN6BUV\nIXmceir06lV/TyyRAlZQ1Vwi+VbvyPN16+CJJ0LyeOEF6NgxJI/HH4eDDtJ4DpEUKZlIy/T227D5\n5uEO5JZbwohxEWk0JRNpmQ47DMaNy3cUIs2G+jJK81bfpIWb6O8okUxSMpHm6bvv4NJLQ88sEck6\n/XkmzU9ZGfziF9CrF8X9+hFLsIaD5mkSySx1DZbmY+VKGDwYnn4ahg+Hk07Kd0QiTVamuwarmkua\nh9degwMPhOXLYcYMJRKRHFM1lxS2Vavg2mvD+JC//jV09RWRnNOdiRSuN98MAwy//DLcjSiRiOSN\n7kyk8KxeDddfD48+CvfdB6efnu+IRFo8JRMpLG+/DRdcEO5Ipk+HbbbJd0TNzoABQykvr9igvEuX\nYkaMGJyHiFJT6PFnU33XJhOUTKQwVFTA738Po0bBvffCmWfmO6KUFNIXXHl5BRMnxhJsSVTW9BR6\n/NlU+9rcmNFzK5lI0zd5Mpx/Puy/f7gb2W67fEeUsvq+4KqqYqxcCWvX1v9Yt67h7Zne74MPEn+G\nqVOhZ8/sXqdM+PDDxOWFEn821XdtMkHJRJquNWsgFoOHHgprqZ91Vr4jSpp7WFfr44/DI8Es+AC8\n/jpsuy384AeJH5tsUv+2VPbbbLPk9120CL75ZsNY99wTbr89u9ctEy69FKZN27C8UOLPpvquTSYo\nmUjTNGVKuBvZe+9wN7L99vmOKCF3+Oyz9Ukj/gGw337QtesGS8LX6NUrDNhvSm69NXH5llsWxl/2\nW25Zf3khxJ9N9V2bTFAykaZlzRq46SZ44AH485/h7LObxNoi7rBgQUgSH31UO2kUF69PGgceGELu\n2jXccVSHXloKc+bk9SOIZJWSiTQdU6fCeefBHnuEe/Edd0z7lKk2fFdVwaefbniXMXMmtGsXkkTX\nrlBSEm6c9t23+XUo69KlmESN1aG86Sv0+LMp/tokWDMuLZqbS/Lv++/DAlX/7//BnXfCz3+esbuR\n0tJYwobvI4+M8cADsQ2SxqxZ0KHD+qRRfcex776w1VaNj6OQenNJy6Ble6WgbbAm+3ffwaxZFLdv\nz+Bp02CnnXISx+uvw/HHr08aRx8NgwbBPvuEO5BMU8KQ5k7JRHKq3jXZDzkkK4mkvhvWww8PCUVE\nMkNzc0nTkIVG9lmz6u8G2apVxt9OpEXLejIxs95mNsvMys3smgb2O9TM1prZT1M9VgrIt99m/S2+\n/x5uvhmOOCL0qBKR7MtqNZeZFQH3AkcDnwFTzOxZd5+VYL+hwEupHisF5PHHw+y+WTR5Mlx4IXTs\nGDqH3XxzMbvsEttgP/XsEcmsbLeZlACz3X0egJmNBk4F6iaEQcATwKGNOFaaOvfQW2vEiDAQ4733\nMv4WK1aEiYQffzwMT/nZz0LNmRq+RXIj28lkZ2BB3OuFhCRRw8x2Ak5z96PMrCSVY6UArFkTbhVm\nzYLJkykeMoTY5ptvsFs6a7K/8AJcfDEcdVSYe6hDh3QCFpHGaAq9ue4C0m4PicViNc9LS0spLS1N\n95SSri+/DGuN7LhjGCHVti2DR4zI2Om/+AJ+/etQtfXAA3DssRk7tUizU1ZWRlkW5+7J6qBFM+sB\nxNy9d/R6MODuPixun7nVT4FtgJXAAOCLjR0bdw4NWmxqPv4YTj4ZzjknTI9SlLm+Hu5hJvqrrw4D\n5mOx+ue+EpHECm3Q4hSgs5l1AhYD/YCz43dw9z2qn5vZw8BYdx9jZq02dqw0US+9BOeeC3fcEX5m\n0Ny5MHAgfP11qN465JCMnl5EGimrXYPdvRK4DBgPfASMdveZZjbQzAYkOmRjx2YzXsmA++4Lk1Y9\n9VRGE8m6dWH68JISOO44eOcdJRKRpkRzc0lmrFsHV14JL78Mzz0XJmvMkPffD234W20Fw4dD584Z\nO7VIi1Vo1VzSEixfDv36QWUlvPVWejMixlm9Gm68MayNNWxYuOFpArPRi0gCmk5F0vPJJ2Giq913\nD40YGUokr7wCBxwQpoOfMQMuuECJRKQp052JNN5bb8EZZ8C118Jll2Xk237pUrjqKpgwITS/nHJK\nBuIUkazTnYk0ziOPwGmnhTqoQYPSTiTuYfT6/vuHbr4ffqhEIlJIdGciqamqCgM7/vGPUBe1//4p\nHZ5okaiKCli8uJjNNx/Mk0/CYYdlMF4RyQklE0ne6tWhFXz+fJg0CbbfPuVTlJdXJFz5sFOnGFOn\nQps26YcpIrmnai5JzpIlUFoaFgJ59dVGJZKG7LabEolIIVMykY2bPh1+9CM46aTQVlLc+OnbKzZc\nBl1EmgFVc0nDnn8+9Mu9554wlqSRFi2CP/4R3n03g7GJSJOhOxNJzD0sDHLRRTBmTKMTyeefh4Hx\nBxwAm24apkMRkeZHdyayobVrw7iRt98Oj06dUj7F11/DbbeF9bD+53/go4/CTPQDBhTTunVsg/21\n8qFIYdPcXFLbsmXQt29oDX/0UWjXLqXDv/kG7rwzDDjs2xeuuw523TVLsYpIo2luLsmYoQMGUFFe\nvr5g9WqYMYPiPfdk8LRpoedWklasgLvvDo+TT4YpUzI616OINHFKJi1YRXk5sYkTNyiPdeiQdCJZ\nuTLchdxxBxxzDLz5JqSxAq+IFCglE2mUigr461/DbL5HHBEGw++3X76jEpF8UTKRlHz/PTz4INxy\nS1ic6sUX4aCD8h2ViOSbkokkZe3asO76H/4A++wTFlJUN18RqaZk0pItX77RXSor4V//CotUdewY\nBsAffngOYhORgqJk0lKtW0fxvHnE9t0Xttuu1qbiLl2oqoInnoAhQ6BDB/jb3+Coo/IUq4g0eUom\nLdV993Hqpu25eNu+wPqu5u7Q5tNiHj0YWrcOg+CPP16rHIpIwzRosSVavBi6daP/7j/jH1Pu3WDz\nZpvF+Ne/YpxyipKISHOV6UGLmpurJbr6arjwQua33Sbh5kMOgT59lEhEJHlKJi3Na6/BxIlw/fX1\n7qIkIiKpqrfNxMymJ3H8l+5+dEM7mFlv4C5C4nrQ3YfV2d4H+ANQBawFrnT3N6NtnwLLq7e5uzqj\npmPtWrj00tAQstlm+Y5GRJqRhhrgWwEnNrDdgDENndzMioB7gaOBz4ApZvasu8+K2+3f7j4m2v8A\n4HFg32hbFVDq7ssa/BSSnHvvhZ12gp/+FIDvvstzPCLSbDSUTAa6+7yGDjazSzZy/hJgdvV5zGw0\ncCpQk0zcfVXc/psTEkjNW6CquMz47LMwbP2tt8CMqipYuLCYvfaKsdNOtXfVdPAikqp6k4m7v2Fm\nrYBR7v7z+vbZyPl3BhbEvV5ISDC1mNlpwJ+AbYGT4t8CmGBmlcAId//bRt5P6nPVVTBwYM0sjKNG\nQceOg3n77ZQmBxYRSajBcSbuXmlmncystbt/n60g3P0Z4BkzOwK4GTg22nS4uy82s20JSWVmfQks\nFovVPC8tLaW0tDRb4RaeV1+FN94IK1UBS5fC4MHw3HNKJCItRVlZGWVlZVk7/0bHmZjZKEIbxhhg\nZXW5u9+50ZOb9QBi7t47ej04HFq7Eb7OMXOAQ919aZ3yIcCKRO+rcSYNWLsWDjwwVHGdfjoAF18M\nRUVh6ngRaZnysTjWnOhRBGyR4vmnAJ3NrBOwGOgHnB2/g5nt6e5zoufdgdbuvtTM2gJF7v6dmW0G\nHAfcmOL7y913h2V3TzsNCItWPfMMfPxxnuMSkWYlmWTylLvPaMzJo2qyy4DxrO8aPNPMBobNPgI4\nw8z6A98Dq4GzosO3B542M4/ifMTdxzcmjhZr4UIYOjSs425GZWW4Kxk2DLbeOt/BiUhzkkw11+tA\nG+DvhC/0jU81m2Oq5qpHv36w115h3njg/vth9OgwZlEDE0VatkxXcyU1N5eZdQEuAPoC7wB/b0p3\nCUomCbz8Mlx4IXz0EbRty+efw/77h7b4/ffPd3Aikm95SSbRG7cCTgPuAb4ljAG51t2fylQwjaVk\nUsf334dG92HDwiRbwHnnhZnmb7stz7GJSJOQ8wZ4M+tGuCs5CZgAnOLuU81sJ+BtIO/JROq46y7Y\nc0845RQgTMf1yiswc2ae4xKRZiuZNpOJwIPA/7n76jrbznX3f2QxvqToziTOggVw8MEweTLsuSdr\n14aXsRiceWa+gxORpiJn1VxmNgJ4kTB31opMvWE2KJnEOess6No1ZA/g9tvh3/+GF19Uo7uIrJfL\nZPIj4ATCJI3fE7r3jnP3DzL15pmiZBIZPx5+9avQ6L7ppixcCAcdFHoG77VXvoMTkaYkX725OhAG\nDZ4AdAOmEhLL45kKJB1KJsCaNdCtG9xxB5x8MgB9+4ablBs11FNE6sjHCHjc/Wvg0eiBmR0C9M5U\nEJIBd94Je+9dk0heegmmTg0TOoqIZFsyDfBbAf2B3YhLPu5+eVYjS0GLvzOZPx+6dw9zpey+OxUV\ncMABYSaVExtakUZEWqx83Jm8AEwCZlB7rRFpKq68Ei6/HHbfHYBbbw3JRIlERHIlmTuTqe7ePUfx\nNEqLvjMZNw4uuww+/BCKi5kzB0pK4P33oWPHfAcnIk1VzhvgzexK4DvgOWBNdXndKeLzqcUmkzVr\nwtwo99wDJ5yAe2gyOfJIuOaafAcnIk1ZPqq5vgduA64jrHxI9HOPTAUhjXT77SGZnHACEKaW/+QT\nePrpPMclIi1OMncmc4ESd/8qNyGlrkXemXz6Kfzwh/Duu7DbbqxcGboBjxwJWmRSRDYm03cmRUns\n819gVabeUDLkyivDY7fdgDDL/I9/rEQiIvmRTDXXSmCamb1K7TaTJtM1uMV54YUwyn30aCCsmvjg\ngzCjUUuYiYikL5lk8kz0kKagogIGDQoLuLdpgztcein8/vewww75Dk5EWqqNJhN3H5mLQCRJt94a\nJtzqHSYg+Ne/YPnysByviEi+NDhrsLsPaPDgJPbJhRbTAD93bhhEMnUqdOzI8uWw777w1FPQo0e+\ngxORQpLLrsGnmVlFQ7EAR2UqEEnCr38Nv/1tzWjEG24I40qUSEQk3xpKJlclcfzrmQpENmLsWPjP\nf+D//g8INyePPRYa30VE8i3pNeCbsmZfzbV6Ney3HwwfDsceS1UV9OwJF10Ev/xlvoMTkUKUj3Em\nkm/DhoUBisceC4RuwEVFcMEFeY5LRCSS9WRiZr3NbJaZlZvZBjNGmVkfM/vAzN43s3fM7PBkj20R\n5syBe+8N65UAX30F110H998fEoqISFOQ1WouMysCyglL/34GTAH6ufusuH3auvuq6PkBwOPuvm8y\nx8ado3lWc1XP3NirF1x9NQAXXgibbw533ZXn2ESkoOW8msvMJkQLZFW/3trMXkry/CXAbHef5+5r\ngdHAqfE7VCeSyOasXzNlo8c2e2PGhO7Av/41AG+9BS++CDfdlOe4RETqSKaiZBt3/6b6hbsvA7ZL\n8vw7AwviXi+Mymoxs9PMbCYwFvhFKsc2W6tWwRVXhCqu1q1Ztw4uuSRMFNyuXb6DExGpLZnpVKrM\nrKO7zwcws06sn4o+I9z9GeAZMzsCuBk4NtVzxGKxmuelpaWUFvqMh3/6UxhAcvTRQJg9pUMH6Ncv\nz3GJSEEqKyujrKwsa+dPZgr63sAIYCJhoOKPgQHuvtGqLjPrAcTcvXf0ejDg7j6sgWPmAIcCXZI9\ntjm0mQwdMICK8vLwYtWqsFTiD39I8QEH0D82gm7d4I03YJ998huniDQPOV8cy93HmVl3oHqc9a9T\nWNtkCtA5uptZDPQDzo7fwcz2dPc50fPuQGt3X2pmGz22OakoLyc2cWLtwkmTiLVpw//+LwwYoEQi\nIk1XMtVcAJXAF0Ax0DXKaK9t7CB3rzSzy4DxhPaZB919ppkNDJt9BHCGmfUnrOi4GjiroWNT/HwF\nb9kyeOtTeOCBfEciIlK/jSYTM7sQuALYBZhGuEN5G/hJMm/g7uOAveuUDY97fitwa7LHtjSzZ8Pd\nj0LbtvmORESkfsn05rqC0IYxz92PAg4Gvmn4EElZPW0+m24KffrkOBYRkRQlk0wq3L0CwMzaRIMG\nW/TdQlYsXpywuHNnsIw1kYmIZEcybSYLo0GLzwATzGwZMC+7YbUwc+ZQNOcTfrb5jqxp1RqAlSuh\nVSsoWr5qIweLiORfMr25To+exqJ14LcExmU1qpakshL696f97kfz+Jw6l3Ud9NohlpewRERSkWxv\nLgDcfeLG95KU3HorFBfz5M4/gjn5DkZEpHE072w+TZsWZgN++GFcDSMiUsCUTPKlogLOPRfuuKNm\nGV4RkUKlZJIvN9wAXbqEhEJoOhERKVQptZlIhrz2GjzyCEyfXtPvd9myYrbbLsa++9betUuX4jwE\nKCKSGiWTXPv2WzjvPBgxArbZBoCyMli2bDAzZ0L79vkNT0SkMbK60mKuFNSswb/8ZVhv929/A2DF\nCujWLSxbctJJeY5NRFqMnM8aLBk0Zky4DZk2rabo6qvhqKOUSESksCmZ5MqXX8KvfgWPPw5bbAHA\nhAnw/PMwY0aeYxMRSZOquXLBHX7609B7a1hY22v5cjjggDC1/HHH5Tk+EWlxVM1ViEaOhLlzYfTo\nmqIrr4QTT1QiEZHmQckk2+bNg6uugpdfhjZtgFC1VVYGH3yQ39BERDJFySSbqqrg/PNDMunWDYCl\nS2HgQPjnP2uaTkRECp5GwGfTXXfBunXw29/WFF1+eWg+KS3NX1giIpmmO5Ns+egj+NOfYPLksDAJ\n8PTT4WVcz2ARkWZBySQbvv8+zLn1pz/BHnsAoWfwJZfAE0/AZpvlOT4RkQxT1+BsuP760Lo+ZkzN\n3FtnnQWdOsFtt+U5NhER1DW46Zs0KQwemTatJpE89lgYmDhyZJ5jExHJEiWTTFq5MlRv3X8/7LAD\nAEuWhEb3sWNh003zHJ+ISJZkvTeXmfU2s1lmVm5m1yTYfo6ZfRA93jCzbnHbPo3K3zezd7Ida9qu\nugp69gzdtQgD33/1K7jwQigpyXNsIiJZlNU7EzMrAu4FjgY+A6aY2bPuPitut7nAke6+3Mx6AyOA\nHtG2KqDU3ZdlM86MGDcujEacPr2m6J//DAPfH3ssj3GJiORAtqu5SoDZ7j4PwMxGA6cCNcnE3SfF\n7T8J2DnutVEIY2GWLg23H6NGwZZbArBoURhe8tJLNQPfRUSarWx/Ue8MLIh7vZDayaKuC4EX4147\nMMHMppjZRVmILzMuuQT69oWf/AQI1VsXXQSXXgoHH5zn2EREcqDJNMCb2VHABcARccWHu/tiM9uW\nkFRmuvsbiY6PxWI1z0tLSynN1RDzRx8NVVsPP1xT9NBD8PnncO21uQlBRGRjysrKKCsry9r5szrO\nxMx6ADF37x29Hgy4uw+rs1834Emgt7vPqedcQ4AV7n5ngm35GWeyaFG49XjxRTjkECDM6/jDH8Kr\nr8L+++c+JBGRZGR6nEm2q7mmAJ3NrJOZtQb6AWPidzCzjoREcm58IjGztma2efR8M+A44MMsx5s8\nd/jFL2DQoJpEUlUVVuX9zW+USESkZclqNZe7V5rZZcB4QuJ60N1nmtnAsNlHADcA7YH7zcyAte5e\nAmwPPG1mHsX5iLuPz2a8DRk6YAAV5eXrCxYtgiVLKN51VwZHRcOHhzXdr7oqLyGKiOSNplNJUqy0\nlNjEiRuW9+pFrKyMuXPhRz+C11+HffbJaigiImkrtGquFqGqCi64AH73OyUSEWmZlEwy4C9/gcpK\nuOKKfEciIpIfTaZrcKFatQruuxnefrtm2RIRkRZHySRZ9bTJzJoFQ/4InTvnOB4RkSZEySRJxUVF\nxDbbLHQDjqaWnz8fllZ04ZJL8hyciEieqTdXMr76Kgwcee65MCIR+Phj6NULpkyB3XbL3luLiGRD\npntzKZklg3UHAAAQ6UlEQVQk4/zzmfDef7mlwzFAqPGaOhV23BFKS4sZMWLwRk4gItK0aKXFXPv3\nv6GsjDt3PYeJE2O1Ns2eDTvtFEt4mIhIS6KuwQ1ZtSqsbnX//axu1Trf0YiINFlKJg35wx9CG8mJ\nJ+Y7EhGRJk3VXPWZPh0efLDWyokiIpKYkkkilZVhdas//hF22AEIU6aISP7ttttuzJs3L99hFIxO\nnTrx6aefZv19lEwSue8+2HTTMJ98pKKimPbtYxxwQO1du3QpznFwIi3bvHnzaA69UHPFLGMdthp+\nn+bwS8lo1+D586F7d3jzTdh7bwCWLAnDTCZN0kh3kXyLurTmO4yCUd/10jiTBDKWTNyhTx8oKYEb\nbqgp/sUvYNttYdiwBo4VkZxQMklNrpKJqrniPfkkzJ0bfkbefRfGjQtzcImISGJKJtW++SbMIf/4\n49A6jClxD0U33wzt2uU5PhGRJkzjTKpdc02o4jr88Jqi0aOhogLOPz9/YYlIy1BUVMTcuXMBuOCC\nC/j973+f54hSo2QCYa3d55+HoUNrilatCvnl7ruhSFdJRLIs1V5XZWVlFBUVcdttt2UpotSommvN\nGhgwAO65B7bcsqb41luhZ0844og8xiYiSRswYCjl5RUblHfpkvxkrJk4R2Ol2qlg1KhRdOjQgVGj\nRnHVVVdlKarkKZkMHRq6AJ9+ek3R/PlhKd73389jXCKSkvLyig0mYw0SlWXvHAsXLuSKK67g9ddf\nx905++yzueeeewB46KGHuP322/n8888pKSlh+PDhdOzYMelzV1u1ahVPPPEEDzzwAP3792fq1Kl0\n79495fNkUsuuwJk5M2SNe++tWfAKQvXWoEHQiN+xiLRgVVVVnHzyyey+++7Mnz+fRYsW0a9fPwCe\nffZZhg4dyjPPPMOXX37Jj3/8Y84+++xGvc+TTz7JFltsQd++fTnuuOMYOXJkJj9Go7TcZFJVFaq3\nYjHYZZea4jfeCOMVr746f6GJSOZMnBj+VkzmMXFieu/1zjvvsHjxYm699VaKi4tp3bo1PXv2BGD4\n8OH87ne/o0uXLhQVFTF48GCmTZvGggULUn6fUaNG0a9fP8yMc845h9GjR1NZWZle8GnKejIxs95m\nNsvMys3smgTbzzGzD6LHG2bWLdlj0/LAA7B2LVx8cU1RVVXoCjxsGLRtm9F3E5E86dUrdPNP5tGr\nV3rvtWDBAjp16kRRgl478+bN44orrqB9+/a0b9+eDh06YGYsWrQopfdYuHAhr776Kueccw4Affr0\nYfXq1Tz//PPpBZ+mrLaZmFkRcC9wNPAZMMXMnnX3+CGAc4Ej3X25mfUGRgA9kjy2cRYvhuuug1de\ngVataor//ncoLoborlREJCW77ror8+fPp6qqaoOE0rFjR66//vpGV21VGzVqFO7OKaecUtNov2bN\nGkaOHEmfPn3SOnc6st0AXwLMdvd5AGY2GjgVqEkI7j4pbv9JwM7JHttoV1wRqrjiZm389tuQX8aO\nrdV8IiIFIky6GqunPDfnKCkpYccdd2Tw4MHEYjFatWrFe++9R8+ePRk4cCA33HADBx54IF27dmX5\n8uVMmDCBM888M+n4ICSTWCzGwIEDa8omT55M3759WbZsGVtvvXVK58uUbCeTnYH4CsGFhCRRnwuB\nFxt5bHLGjg3dtOo0WN1yC5xwQlgLS0QKTya67qZ7jqKiIsaOHcugQYPo2LEjRUVFnHPOOfTs2ZPT\nTjuNlStX0q9fP+bPn8+WW27JscceW5NMkhlnMnnyZObPn88ll1xChw4daspPOeUU9tprLx599FEu\nueSStD5DY2V1okczOwM43t0HRK//Byhx98sT7HsUoVrrCHdfluKxPmTIkJrXpaWllJaWbhjQihWw\n336hPusnP6kp/u9/oUcP+PDDmuVLRKSJ0kSPqam+XmVlZZSVldWU33jjjYUza7CZ9QBi7t47ej0Y\ncHcfVme/bsCTQG93n5PKsdG2hLMGDx0wgIry8vUFs2dDZSXFffoweMSImuLTToPDDgtdgkWkaVMy\nSU1zmTV4CtDZzDoBi4F+QK3WJzPrSEgk51YnkmSP3ZiK8nJiCfr6xeISzIQJMGMGPPZYKmcWEZF4\nWU0m7l5pZpcB4wndkB9095lmNjBs9hHADUB74H4LlYZr3b2kvmMzGd+6dXDllXDHHdCmTSbPLCLS\nsmR9OhV3HwfsXadseNzzi4CLkj02k4YPh+23h1NPzdY7iIi0DC12bq6vv4Ybb4SXX1ZXYBGRdLXY\n6VRiMejbt9ZQExERaaRmfWdS3KVLwrk+V23Thcceg48/znVEIiLNU1a7BudKfV2DE3GH44+Hk0+G\nyzcYsSIiTZ26BqcmV12DW1w113PPwcKFteZ3FBHJOy3bW0DWrIHf/Ab+/Gf4wQ/yHY2IyHrJLts7\ncuRINtlkE9q1a0e7du3o3Lkzf/3rX7Mc3cY16zaTuu65B/bZJ1RziUjzssGMF5HiLl1qzXiR7XM0\nVipVdz179uS1114DYNq0aRx55JEcdthhHHjggdkKb6OadTKJX8/5++9hyhQ4+GAYMCD76zmLSG7V\nO+NFjs+Ri2V74x100EHsu+++zJw5U8kkWxKt5zxlCrRtG0u4v4hIOqqX7T3mmGN45JFHKCoq4t13\n3wXWL9v73HPP0blzZ4YOHcrZZ5/Nm2++mdZ7TpkyhdmzZ/PDPE953qLaTESkBcrhur25Wrb37bff\npn379rRr144ePXpw7rnn0rlz57RiT5eSiYg0bzlctzcXy/YCHHbYYSxdupRvv/2WJUuW8OGHH3Ld\nddelFXu6lExERDIkftneujp27Mjw4cNZunQpS5cuZdmyZXz33Xf06NEjrffcdtttOeOMMxg7dmxa\n50lXs24zEZGWo74ZL4q7dMnZOXKxbC/U7vn19ddf8/TTT7P//vunfJ5MatbJJBNrQotIYchE1910\nz5HtZXurTZo0iXbt2gHQtm1bjjnmGO666660Yk9Xi5tORUQKm6ZTSY2mUxERkYKhZCIiImlTMhER\nkbQpmYiISNqUTEREJG1KJiIikrZmPc5ERJqfTp06pTQmo6Xr1KlTTt4n6+NMzKw3cBfhLuhBdx9W\nZ/vewMNAd+Bad78zbtunwHKgCljr7iX1vIfGmYiIpKCgxpmYWRFwL3A8sB9wtpntU2e3r4FBwG0J\nTlEFlLr7wfUlkkJXVlaW7xDSovjzS/HnV6HHn0nZbjMpAWa7+zx3XwuMBk6N38Hdv3L394B1CY63\nHMSYV4X+j1Hx55fiz69Cjz+Tsv1FvTMQP1n/wqgsWQ5MMLMpZnZRRiMTEZGMaeoN8Ie7+2Iz25aQ\nVGa6+xv5DkpERGrLagO8mfUAYu7eO3o9GPC6jfDRtiHAivgG+GS3m5la30VEUpTJBvhs35lMATqb\nWSdgMdAPOLuB/Ws+mJm1BYrc/Tsz2ww4Drgx0UGZvCAiIpK6rCYTd680s8uA8azvGjzTzAaGzT7C\nzLYH3gW2AKrM7AqgK7At8HR017EJ8Ii7j89mvCIi0jjNYj0TERHJr4Ludmtmvc1slpmVm9k1+Y6n\nPmb2qZl9YGbvm9k7UdnWZjbezP5jZi+Z2ZZx+//OzGab2UwzOy4P8T5oZp+b2fS4spTjNbPuZjY9\n+v3kZBm4emIfYmYLzWxq9OjdFGOP3ncXM3vFzD4ysxlmdnlUXijXv278g6LygvgdmFkbM5sc/V+d\nEbXVFsT1byD23Fx7dy/IByER/hfoBPwAmAbsk++46ol1LrB1nbJhwNXR82uAodHzrsD7hKq93aLP\naDmO9wjgIGB6OvECk4FDo+cvAMfnKfYhwG8S7LtvU4o9eq8dgIOi55sD/wH2KaDrX1/8hfQ7aBv9\nbAVMIoyXK5Trnyj2nFz7Qr4z2eiAyCYk0eDLU4GR0fORwGnR8z7AaHdf5+6fArMJnzVnPHS/Xlan\nOKV4zWwHYAt3nxLtNyrumKypJ3aI69wR51SaUOwA7r7E3adFz78DZgK7UDjXP1H81WPLCuV3sCp6\n2obwResUzvVPFDvk4NoXcjJJd0BkLsUPvrwwKtve3T+H8B8Q2C4qr/u5FtE0Ptd2Kca7M+F3Ui3f\nv5/LzGyamT0QV0XRpGM3s90Id1mTSP3fS94/Q1z8k6OigvgdmFmRmb0PLAEmRF+qBXH964kdcnDt\nCzmZFJLD3b07cCJwqZn9mPV/MVQrtJ4QhRTv/cAe7n4Q4T/ZHXmOZ6PMbHPgCeCK6C/8gvr3kiD+\ngvkduHuVux9MuCMsMbP9KJDrnyD2ruTo2hdyMlkEdIx7vUtU1uS4++Lo55fAM4Rqq88tdIsmuq38\nItp9EbBr3OFN5XOlGm+T+Rzu/qVHlb/A31hfbdgkYzezTQhfxP9w92ej4oK5/oniL7TfAYC7fwuU\nAb0poOsPtWPP1bUv5GRSMyDSzFoTBkSOyXNMGzCzttFfadj6wZczCLGeH+12HlD9pTEG6Gdmrc1s\nd6Az8E5Ogw6M2vWsKcUbVQUsN7MSMzOgf9wx2VYr9ug/f7WfAh9Gz5ti7AAPAR+7+91xZYV0/TeI\nv1B+B2a2TXU1kJltChxLaPdp8te/nthn5ezaZ7t3QTYfhL8Y/kNoOBqc73jqiXF3Qk+z9wlJZHBU\n3h74dxT/eGCruGN+R+hZMRM4Lg8x/wv4DFgDzAcuALZONV7gkOgzzwbuzmPso4Dp0e/hGUL9d5OL\nPXrfw4HKuH8zU6N/5yn/e8nT9a8v/oL4HQAHRDFPi+K9Lipv8te/gdhzcu01aFFERNJWyNVcIiLS\nRCiZiIhI2pRMREQkbUomIiKSNiUTERFJm5KJiIikTclERETSpmQizZaZVUbrN7wf/bw6Kn/VzLpn\n+L2uMLPiuNfPmVm7DJz3YTOba2YDUjzuFTNbkenPKVKfbK8BL5JPKz1MsJkLvwb+AVQAuPvJGTz3\n/7r7U6kc4O4/MbNXMhiDSIN0ZyLNWaI1HGrvYHasmb1lZu+a2WPRXGrHm9njcfv0MrMx0fP7zeyd\nOivZDQJ2Al41s5ejsk/MrH30/DfR/tPN7IqorJOZfWxmI8zsQzMbZ2Ztkoj34SiGt83sv1FsD0bn\neijVzy+SKUom0pxtWqeaq2/8RjPrAFwPHO3uPwTeA35DmIOpJJosD+BnhMXXAK519xLgQKDUzPZ3\n978QZlUtdfejo/08eo/uhIkBDwUOAy4yswOjfToDf3H3/YHlwBlJfq6t3P2wKNYxwB3u3hXoZmbd\nkjyHSEapmkuas1UbqebqQVh29c1odtQfAG+5e6WZjQNOMbMngZOAq6Jj+pnZRYT/OztEx3/IhrMs\nVzsCeNrdKwDM7Cngx8BY4BN3nxHt9x5h6dRkjI1+zgCWuPvH0euPonNMT/I8IhmjZCItmQHj3f3n\nCbY9BlxGWAJ4iruvtLBy4G+BQ9z9WzN7GChOcGyy1sQ9r0zhXNXHVdU5RxX6Py15omouac421mYw\nCTjczPaEmrVn9oq2TQS6AxexvoqrHfAdsCJaKOmEuHN9G22v+96vA6eZWXG0ns3pUVky8SVD7SLS\nJOivGGnOis1sKuEL14Fx7n5t9Bx3/8rMzgcejRq/ndCGMtvdq8zsOUJ7R/9o/+lmNo2w9sMC4I24\n9/obMM7MFkXtJtXv8b6Z/Z2wmJsDI9z9AzPrROOWfm1o+VitJyF5o/VMRJqwqCrtOXd/shHHvgr8\n1t2nZj4ykdpUzSXStC0HbmrMoEXCKp9rsxKVSB26MxERkbTpzkRERNKmZCIiImlTMhERkbQpmYiI\nSNqUTEREJG3/HynBatZzJE2eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fabb1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(elev_mec,ac['cellA'],'bs-',label='cell A')\n",
    "plt.plot(elev_mec,ac['cellB'],'rs-',label='cell B')\n",
    "plt.xlabel('Elevation [m]')\n",
    "plt.ylabel('ac [m/yr]')\n",
    "plt.title('Accumulation')\n",
    "plt.legend(loc='lower right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SMB at point P is now estimated by interpolating the regression parameters for ablation, and linear interpolation of accumulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a_P = 0.5*(intercept['cellA']+intercept['cellB'])\n",
    "b_P = 0.5*(gradient['cellA']+gradient['cellB'])\n",
    "ab_P = a_P + b_P * elev_P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_interpA = interp1d(elev_mec, ac['cellA'])\n",
    "yA = y_interpA(elev_P)\n",
    "y_interpB = interp1d(elev_mec, ac['cellB'])\n",
    "yB = y_interpB(elev_P)\n",
    "ac_P=0.5*(yA+yB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'method1': -2.5628565450650003,\n",
       " 'method2': -2.5628565450649998,\n",
       " 'method3': -2.5680887639600036}"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smb_P['method3'] = ac_P - ab_P\n",
    "smb_P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}