gustavo-marques/vgrid_65L_6000m.ipynb

## vgrid_65L_6000m.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vertical grid CESM/MOM6 \n",
    "\n",
    "### Requirements:\n",
    "* NK = 63\n",
    "* MAXIMUM DEPTH = 6000 m\n",
    "* MINIMUM DEPTH = 10 m\n",
    "* The sum of the first four levels must be equal to the MINIMUM DEPTH\n",
    "* $\\sum{dz}$ = MAXIMUM DEPTH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First, let's explore the 'OLD' vertical grid used in CESM/MOM6 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "old_vgrid = xr.open_dataset('/glade/p/cesmdata/cseg/inputdata/ocn/mom/tx0.66v1/ocean_vgrid_190320.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.Dataset&gt;\n",
       "Dimensions:  (nk: 63)\n",
       "Dimensions without coordinates: nk\n",
       "Data variables:\n",
       "    dz       (nk) float32 ...\n",
       "Attributes:\n",
       "    title:    Vertical resolution for the MOM6 tx0.66 grid using Z*\n",
       "    Author:   Gustavo Marques (gmarques@ucar.edu)\n",
       "    Created:  2019-03-20T09:38:19.474327</pre>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (nk: 63)\n",
       "Dimensions without coordinates: nk\n",
       "Data variables:\n",
       "    dz       (nk) float32 ...\n",
       "Attributes:\n",
       "    title:    Vertical resolution for the MOM6 tx0.66 grid using Z*\n",
       "    Author:   Gustavo Marques (gmarques@ucar.edu)\n",
       "    Created:  2019-03-20T09:38:19.474327"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "old_vgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.DataArray &#x27;dz&#x27; (nk: 63)&gt;\n",
       "array([  2.5 ,   2.5 ,   2.5 ,   2.5 ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.25,  10.5 ,  11.  ,  11.75,  12.5 ,  13.25,  14.25,  15.5 ,\n",
       "        17.  ,  18.75,  20.5 ,  23.  ,  25.5 ,  28.75,  32.5 ,  36.75,  42.  ,\n",
       "        48.  ,  55.25,  63.75,  73.75,  85.25,  98.5 , 113.25, 129.75, 147.  ,\n",
       "       164.75, 182.  , 198.  , 211.75, 223.25, 231.75, 238.25, 242.5 , 245.5 ,\n",
       "       247.25, 248.5 , 249.  , 249.5 , 249.75, 250.  , 250.  , 250.  , 250.  ],\n",
       "      dtype=float32)\n",
       "Dimensions without coordinates: nk\n",
       "Attributes:\n",
       "    long_name:  Nominal thickness of level\n",
       "    units:      dz</pre>"
      ],
      "text/plain": [
       "<xarray.DataArray 'dz' (nk: 63)>\n",
       "array([  2.5 ,   2.5 ,   2.5 ,   2.5 ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.25,  10.5 ,  11.  ,  11.75,  12.5 ,  13.25,  14.25,  15.5 ,\n",
       "        17.  ,  18.75,  20.5 ,  23.  ,  25.5 ,  28.75,  32.5 ,  36.75,  42.  ,\n",
       "        48.  ,  55.25,  63.75,  73.75,  85.25,  98.5 , 113.25, 129.75, 147.  ,\n",
       "       164.75, 182.  , 198.  , 211.75, 223.25, 231.75, 238.25, 242.5 , 245.5 ,\n",
       "       247.25, 248.5 , 249.  , 249.5 , 249.75, 250.  , 250.  , 250.  , 250.  ],\n",
       "      dtype=float32)\n",
       "Dimensions without coordinates: nk\n",
       "Attributes:\n",
       "    long_name:  Nominal thickness of level\n",
       "    units:      dz"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "old_vgrid.dz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2b4b28da1550>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hc9Zn28e8jybIs915kG7kbG+OCMMWQmF5CDISyIQk9IZuLUBbYBVLJZrMvSYBUNqH3ZsBgk1BDMc0UueCCsTGukmxLbrIsq+t5/5hjECBLY6GZM+X+XNdcM+fMGfk+ZvCjc37N3B0REZF9lRF2ABERSU4qICIi0iYqICIi0iYqICIi0iYqICIi0iZZYQeIpz59+nh+fn7YMUREksr8+fO3uHvfL+5PqwKSn59PYWFh2DFERJKKma1rbr9uYYmISJuogIiISJuogIiISJuogIiISJuogIiISJskTAExsyFm9qqZLTezZWZ2RbD/BjMrNrNFwePkJp+53sxWmdkKMzshvPQiIuknkbrx1gNXu/sCM+sKzDezl4L3/uDuNzU92MzGAd8GxgODgH+Z2Wh3b4hrahGRNJUwBcTdNwIbg9cVZrYcyGvhI6cCj7p7DbDGzFYBU4F5MQ8rItIOGhqd7btrqaypp6qugd21DVTXRp6r6hqoCp73bPMVlt+46Ihh9MjNbsf0CVRAmjKzfGAy8C4wDfixmZ0HFBK5StlOpLi80+RjRTRTcMzsEuASgKFDh8Y0t4gIQHVdAx9u3MmGbbvZuquWbZW1bK2sZVtlzaevt1fWsqOqbp9qglnbM5150JDULyBm1gV4ErjS3Xea2d+AXwMePN8MXAQ091f5pf8U7n47cDtAQUGBVs8SkXbl7nxcuov31mxjSVE5i4vLWbm5gobGz/65yTDo1Tn708f+A7rRs3MHenXuSO/O2XTNyaJTh0w6ZWd++pybnUlOh0xysyPv5XTIwL5KBYmBhCogZtaBSPF4yN1nAbj75ibv3wH8I9gsAoY0+fhgoCROUUUkjVVU1/HWqq3MXVnK3BVllJRXA9AztwMTBvfgmLH9mDC4OyP6dqZ3545079SBjIzE+se/PSRMAbFIab0LWO7utzTZPzBoHwE4HVgavJ4DPGxmtxBpRB8FvBfHyCKSZj7atJO/vrKK55duor7R6doxi2kj+3D5MX2ZNrIPg3t2SrirhFhKmAJCpK3jXGCJmS0K9v0EOMfMJhG5PbUW+CGAuy8zs5nAh0R6cF2qHlgiEgtLisr5yysf8+KHm+nSMYvzD8/nuHH9OWi/nnTITJjREHGXMAXE3d+k+XaNZ1v4zG+A38QslIiktY83V/C/zy7n1RVldMvJ4opjRnHhtPx2b4xOVglTQEREEkVjo3P3W2v43QsryM3O5D9PGMO5h+1Ht5wOYUdLKCogIiJNFG3fzTWPf8A7q7dx7P79ufGMCfTp0jHsWAlJBUREhEh33FkLirlhzjIa3fndmQdy1kGD06pRfF+pgIhI2nN3fvv8Cv4+9xOmDuvFzWdNZEiv3LBjJTwVEBFJe//32if8fe4nfPeQofz3qQeQmYJjNmJBBURE0tq9b63h9y+s4FuT8/j1qQek5IC/WEnfDswikvYeL9zADc98yAnj+/O7Mw9U8dhHKiAikpaeXbKRa59czJGj+vDncyaTlcYDAttKf2MiknbeW7ONKx5dyJShPbnt3IPomJUZdqSkpAIiImmlfHcdVz66kLwenbjrgoPJzVZTcFvpb05E0oa7c92sxZRW1PDkjw6neyeNLP8qdAUiImnj0fc38NzSTVxzwhgmDukRdpykpwIiImlhVWkFv3pmGUeM7MMlRw4PO05KUAERkZRXXdfAZY8sIjc7i1vOnqjuuu1EbSAikvJ++/xHLN+4k7svKKBft5yw46QMXYGISEorXLuNe95aywWH53P02P5hx0kpKiAikrLcnRuf+4h+XTty7Yljw46TclRARCRlvby8lMJ127ni2FF0ytZgwfamAiIiKamh0fndCx8xrE9nzi4YEnaclKQCIiIp6amFxazcvItrjh9DB81zFRP6WxWRlFNd18AfXlrJhLzunDxhQNhxUpYKiIiknAffWUfxjiquPXGslqSNIRUQEUkpO6vruPXVVRwxsg9HjOoTdpyUpgIiIinlztdXs313nbrtxoEKiIikjG2Vtdz55hq+ceBAJgzuHnaclKcCIiIp48F31rG7toErjhkVdpS0oAIiIimhuq6B++etZfqYvozu3zXsOGlBBUREUsKcRSVs2VXL94/QVO3xogIiIknP3bnzzdWMHdCVaSN7hx0nbaiAiEjSe+PjLazcvIvvHzlc4z7iaK/rgZjZM4Dv7X13nxGTRCIi++iON1bTr2tHZkwcFHaUtNLSglI3xS0FYGZDgPuBAUAjcLu7/8nMegGPAfnAWuBsd99ukV8z/gScDOwGLnD3BfHMLCLhW7Gpgjc+3sJ/njCG7CzdVImnvRYQd5+757WZdQKGuvuKGGapB6529wVm1hWYb2YvARcAL7v7jWZ2HXAdcC1wEjAqeBwC/C14FpE0cucbq8npkMF3pg4NO0raabVcm9k3gUXA88H2JDOb095B3H3jnisId68AlgN5wKnAfcFh9wGnBa9PBe73iHeAHmY2sL1ziUjiKq2oZvaiEs46aAg9O2eHHSftRHO9dwMwFdgB4O6LiNxOihkzywcmA+8C/d19Y/BnbwT6BYflARuafKwo2PfFn3WJmRWaWWFZWVksY4tInD0wbx11jY1cdMSwsKOkpWgKSL27l8c8ScDMugBPAle6+86WDm1m35ca/d39dncvcPeCvn37tldMEQlZdV0DD76zjmP378+wPp3DjpOWoikgS83sO0CmmY0ys78Ab8cijJl1IFI8HnL3WcHuzXtuTQXPpcH+IqDpMmODgZJY5BKRxPPSh5vZvruO8w7bL+woaSuaAnIZMB6oAR4GyoEr2ztI0KvqLmC5u9/S5K05wPnB6/OB2U32n2cRhwLle251iUjqm1m4gbwenZg2QlO2h6Wlbrx7jHH3nwI/jXGWacC5wBIzWxTs+wlwIzDTzC4G1gNnBe89S6QL7yoi3XgvjHE+EUkQRdt38+aqLVx+9CgyMjRwMCzRFJBbgltHjwOPuvuyWARx9zdpvl0D4Jhmjnfg0lhkEZHE9uT8YgDOPGhwyEnSW6u3sNz9KGA6UAbcbmZLzOxnsQ4mItKcxkbn8fkbmDaiD0N65YYdJ61FNWzT3Te5+5+BfycyJuQXMU0lIrIX81ZvpWh7FWcV6OojbNEMJNzfzG4ws2XAX4n0wNJ/OREJxczCDXTLyeKE8QPCjpL2omkDuQd4BDjO3dVNVkRCU767jueWbuLbBw8hp0Nm2HHSXqsFxN0P3TMXVhzyiIjs1ZwPiqmtb+TsgiGtHywxlzBzYYmItGZmYRHjBnbjgLzuYUcREnQuLBGRL/qwZCdLiss5W43nCSPh5sISEWnOzMINZGdmcOqkL82ZKiGJphH9c3NhAZcTo7mwRESaU1vfyOxFxRw3vr+mbU8g+zoX1iPATmIwF5aIyN68vrKM7bvrOGOKrj4SSTS9sHYTmQcr1nNhiYg06+lFxfTqnM2Ro7QkQyLZawExs2doZn2NPdx9RkwSiYg0UVFdx0sfbubsgiF0yNSa54mkpSuQm+KWQkRkL15Ytpma+kZOmzwo7CjyBXstIO4+N55BRESaM3tRMUN6dWLK0J5hR5Ev0PWgiCSs0opq3lq1hdMm5RFZc04SiQqIiCSsZz7YSKOjsR8Jaq8FxMweCJ6viF8cEZHPPL2wmAPyujGyX5ewo0gzWroCOcjM9gMuMrOeZtar6SNeAUUkPX1StoslxeWcpquPhNVSL6y/E5lAcTgwn88vN+vBfhGRmJi9sBgz+OZE9b5KVHu9AnH3P7v7/sDd7j7c3Yc1eah4iEjMuDtPLyrh8BG96d8tJ+w4shfRjET/kZlNBI4Mdr3u7otjG0tE0tmC9TtYv203lx09Muwo0oJo1gO5HHgI6Bc8HjKzy2IdTETS1+xFxXTMyuDEA7RsbSKLZjbe7wOHuHslgJn9FpgH/CWWwUQkPdU3NPLPxRs5Zv9+dM3pEHYcaUE040AMaGiy3cDnG9RFRNrN259sZWtlLTMmqvdVoovmCuQe4F0zeyrYPg24K3aRRCSdzfmghK4ds5g+RjPvJrpoGtFvMbPXgCOIXHlc6O4LYx1MRNJPdV0DLyzdxPHjB5DTITPsONKKaK5AcPcFwIIYZxGRNDd3ZRkVNfXMmKSxH8lAc2GJSMKY80EJvTpnc/iI3mFHkSiogIhIQqisqefl5Zs5ecIALRyVJKIZB9LZzDKC16PNbIaZqW+diLSrfy3fTHVdo3pfJZFoyvzrQI6Z5QEvAxcC98YylIiknzmLShjYPYeC/bRwVLKIahyIu+8GvgX8xd1PB8a1dxAzu9vMSs1saZN9N5hZsZktCh4nN3nvejNbZWYrzOyE9s4jIvGzY3ctr39cxikHDiQjQ8PMkkVUBcTMDgO+C/wz2BdV7619dC9wYjP7/+Duk4LHs0GgccC3gfHBZ/7PzNTnTyRJPb90E3UNrttXSSaaAnIlcD3wlLsvM7PhwKvtHcTdXwe2RXn4qcCj7l7j7muAVcDU9s4kIvHxzOIS8nvnckBet7CjyD5otYC4+1x3n+Huvw0a07e4++VxyLbHj81scXCLa8/N0TxgQ5NjioJ9X2Jml5hZoZkVlpWVxTqriOyj0opq5n2ylRkTB2nd8yQTTS+sh82sm5l1Bj4EVpjZf8Y+GgB/A0YAk4CNwM17YjVzrDf3A9z9dncvcPeCvn01NYJIovnn4si651o4KvlEcwtrnLvvJDIH1rPAUODcmKYKuPtmd29w90bgDj67TVUEDGly6GCgJB6ZRKR9zfmghLEDujKqf9ewo8g+iqaAdAjGfZwGzHb3Ovby2357M7OBTTZPB/b00JoDfNvMOprZMGAU8F48MolI+1m3tZKF63dw2mQ1niejaHpT3QasBT4AXjez/YCd7R3EzB4BpgN9zKwI+CUw3cwmESlYa4EfAgSN+TOJ3FKrBy5194bmfq6IJK7Zi0owgxm6fZWUzH3fLybMLMvd62OQJ6YKCgq8sLAw7BgiQmTd82Nunkvfrh157IeHhR1HWmBm89294Iv7o2lE729md5nZc8H2OOD8GGQUkTSypLic1VsqOV23r5JWNG0g9wIvAHuuMVcSGRsiItJmTy8sITszg5MmDGz9YElI0RSQPu4+E2gECG5dqb1BRNqsodF5ZnEJR43tS/dOmps1WUVTQCrNrDdBzyszOxQoj2kqEUlpb3+yhbKKGk6bpNtXySyaXlhXEek2O8LM3gL6AmfGNJWIpLSnFhbTNSeLo8b2CzuKfAXRrIm+wMy+DowhMgJ8RTAWRERkn1XVRtY9P+XAQVr3PMlFO6vuVCA/OH6KmeHu98cslYikrH8t30xlbQOnTtbYj2TXagExsweIzEe1iM8azx1otoCY2Z+j+HN3uvvPog0pIqlj9qJiBnTL4dBhWvc82UVzBVJAZD6saEccngr8opVjrgNUQETSzLbKWl5bUcbFRwzTwlEpIJoCshQYQGQ23Gj8wd3va+mAJtOyi0ga+eeSjdQ3Oqeq91VKiKaA9AE+NLP3gJo9O919RnMHu/sfAcxsiLs3XbMDMxvg7pv2HCMi6eWpBUWM7t+F/Qdq5t1UEE0BuaGNP3uNmT0OXBysqQ6R6eCntPHniUgSW122iwXrd3D9SWO1cFSKiGYgYW6wKuGnD2D/KD63BHgDeMPMRgT79K0RSVNPLigiw9DcVykkmgLyczM7es+GmV1LpKG8Ne7u/wdcDjxjZt8kTuuIiEhiaWx0nlpQzNdG96Vft5yw40g7ieYW1gzgH8EyticCY4N9rTEAd3/LzI4BHgs+KyJpZt7qrZSUV3P9ydHcvJBkEc1I9C1mNgP4FzAfODPKLr0nN/kZG4OrmMPbnFREktYT84vompPFceP6hx1F2tFeC4iZVRC55WTBczYwHDjTzNzdu+3lc1c1ed3cIa9/lcAiklx21dTz/NJNnD4lT1OXpJi9FhB3b2s/uz2fGwMcTGQiRoBvouIhknaeXbKRqroGzpgyOOwo0s6imcrkdOAVdy8PtnsA09396eaOd/dfBce9CExx94pg+wbg8XbKLSJJ4on5RQzv05kpQ3uEHUXaWTS9sH65p3gAuPsO4JdRfG4oUNtku5bIhIwikibWb93Ne2u2ccZBgzX2IwVF0wuruSITzeceAN4zs6eItKGcDrQ4xYmIpJYnFxRhGvuRsqIpBIVmdgtwK5FCcBmR3lgtcvffmNnzwBHBrgvdfWGbk4pIUmlsdGYtLGLaiD4M6tEp7DgSA9HcwrqMyO2nx4i0YVQDl+7tYDNbsOe1u8939z8Fj4XNHSMiqen9tdvYsK2KMw7S1UeqimYcSCWR6dejtb+ZLW7hfQO678PPE5Ek9FjhBrp0zOKE8QPCjiIx0tI4kD+6+5Vm9gzNTEGyt9l4iW60eUPrh4hIstpWWcs/Fm/k3wqGkJsd7cKnkmxa+i/7QPB80778QHdf1/Y4IpIKHnt/A7X1jZx72H5hR5EYamkg4fzgeW784ohIsmtodB56dx2HDu/F6P5a9yOVtdqIbmbTzOwlM1tpZqvNbI2ZrY5HOBFJPq+tKKVoexXnHZYfdhSJsWhuTt4F/AeRrrtquxCRFt0/bx39u3XUxIlpIJoCUu7uz8U8iYgkvXVbK5m7sowrjx1Fh8xoRglIMmupF9aepWdfNbPfA7P4/JroGsshIp/z4DvryMowzpk6NOwoEgctXYHc/IXtgiavHTiadmRmdwOnAKXufkCwrxeRAYz5wFrgbHffbpFJdf5EZM2R3cAFKmgi4aqqbWBmYREnjB9Af606mBZa6oV1FICZDXf3zzWam9nwGGS5F/grcH+TfdcBL7v7jWZ2XbB9LXASMCp4HAL8LXgWkZA8s7iE8qo6dd1NI9HcpHyimX3tPi27u78ObPvC7lP5bALG+4DTmuy/3yPeAXqY2cD2ziQi0XF3Hpi3jtH9u3DIsF5hx5E4aakNZCwwHuhuZt9q8lY3IF7Xp/3dfSN8uixuv2B/HrChyXFFwb6NX/wBZnYJcAnA0KG6LysSC4s27GBJcTm/PnW8pm1PIy21gYwh0ibRg8hqgntUAD+IZagoNPcNbXaddne/HbgdoKCgIJq13EVkH9315hq6dszidK06mFZaagOZDcw2s8PcfV4cMzW12cwGBlcfA4HSYH8RMKTJcYOBkrinExHWba3k2SUb+cHXhtOlo+a9SiettoGEWDwgsp76+cHr84HZTfafZxGHEhmr8qXbVyISe3e8sZqsjAwunjYs7CgSZwnz64KZPQJMB/qYWRGRZXNvBGaa2cXAeuCs4PBniXThXUWkG++FcQ8sImzZVcPjhUV8a0oe/dR1N+0kTAFx93P28tYxzRzrtLColYjEx31vr6W2oZEffC0WPfsl0bXUC+uqlj7o7re0fxwRSRaVNfXcP28dx4/rz4i+XcKOIyFo6QpE8zCLyF498t56yqvq+Pevjwg7ioSkpV5Yv4pnEBFJHrX1jdz15hoOGdaLyUN7hh1HQtJqG4iZ5QAXExlU+GkrmbtfFMNcIpLAnvmghI3l1fzv6RPCjiIhimYqkweAAcAJwFwiYy4qYhlKRBJXY6Nz2+ufMHZAV6aP6Rt2HAlRNAVkpLv/HKh09/uAbwD6tUMkTb38USkrN+/ih18frmlL0lw0BaQueN5hZgcA3YlMry4iaaah0bnphRXk987llAMHhR1HQhZNAbndzHoCPycyAvxD4HcxTSUiCenphcWs2FzB1ceP0YqD0nojurvfGbycC2i0kEiaqqlv4JaXVnJAXje+MUGrJ0h0vbA6AmcQuW316fHu/t+xiyUiiebBd9ZTvKOKG8+YQEaG2j4kuqlMZgPlwHyarIkuIumjorqOW19dxbSRvTlylHpeSUQ0BWSwu58Y8yQikrDueH012yprufbEsWFHkQQSTSvY22ambrsiaaqsooY731zDNyYM5MDBPcKOIwkkmiuQI4ALzGwNkVtYRmRC3ANjmkxEEsJfX/mYmvpGrj5+dNhRJMFEU0BOinkKEUlIa7dU8vB76zm7YAjDNeOufEFL07l3c/edaNoSkbTk7vxyzjKyMzO48thRYceRBNTSFcjDwClEel85kVtXezgaEyKS0p5buom5K8v4+Snj6K/VBqUZLU3nfkrwrIWORdLMrpp6fvXMMsYN7Mb5h+0XdhxJUFEtaWtmB/LlgYSzYpRJREJ2y4srKa2o4e/fO4gsTVkiexHNSPS7gQOBZUBjsNsBFRCRFLSspJx7317DOVOHarEoaVE0VyCHuvu4mCcRkdA1Njo/fWopPXOzufYEDRqUlkVzbTrPzFRARNLAI++vZ9GGHfz0G/vTPbdD2HEkwUVzBXIfkSKyCQ0kFElZpRXV/Pa5jzh0eC9On5wXdhxJAtEUkLuBc4ElfNYGIiIppLHRuXrmB9Q2NPI/p03QSoMSlWgKyHp3nxPzJCISmjvfXM0bH2/hN6cfwMh+GnEu0YmmgHxkZg8Dz9BkOnd14xVJDUuKyvn9Cys4cfwAvjN1aNhxJIlEU0A6ESkcxzfZp268Iimgsqaeyx9dSJ8uHbnxDN26kn0TzZK2F8YjiIjE3y/nLGPt1koe+cGh9MjNDjuOJJlWu/Ga2WAze8rMSs1ss5k9aWaD4xFORGJnzgclPDG/iB8fNZJDh/cOO44koWjGgdwDzAEGAXlE2kLuiWUoEYmtVaW7+OmsJUwZ2oMrjtFMu9I20RSQvu5+j7vXB497AS2KLJKktlXWctG979OxQwZ/Pmey5rqSNovmm7PFzL5nZpnB43vA1lgHa8rM1prZEjNbZGaFwb5eZvaSmX0cPGvSHpFW1NQ38MMHCtm0s5rbzi1gcM/csCNJEoumgFwEnA1sAjYCZwb74u0od5/k7gXB9nXAy+4+Cng52BaRvXB3rn9yCe+v3c7NZ03koP30O5d8NdH0wloPzIhDln11KjA9eH0f8BpwbVhhRBLdra+uYtbCYq46bjTfnDgo7DiSAlpa0vYXLXzO3f3XMciz1z8PeNHMHLjN3W8H+rv7xiDMRjPr19wHzewS4BKAoUM1SErS0z8Wl3DTiys5fXIelx09Muw4kiJaugKpbGZfZ+BioDcQzwIyzd1LgiLxkpl9FO0Hg2JzO0BBQYHHKqBIonp1RSlXzfyAgv16arCgtKuWlrS9ec9rM+sKXAFcCDwK3Ly3z8WCu5cEz6Vm9hQwFdhsZgODq4+BQGk8M4kkg5c+3MylDy1g9IAu3HFeAR2zMsOOJCmkxUb0oKfT/wCLiRSbKe5+rbvH7R9rM+scFDDMrDORKVWWEhmbcn5w2PnA7HhlEkkGzy3ZyI8enM/+g7rx0PcPpWdnjTSX9tVSG8jvgW8Ruf0zwd13xS3V5/UHngouu7OAh939eTN7H5hpZhcD64GzQsonknDmfFDCfzy2iElDenDvhQfTNUeLQ0n7M/fmmwXMrJHIJIr1RBqxP32LSCN6t9jHa18FBQVeWFgYdgyRmJq1oIhrHv+Agvxe3HPBwXTuGM2cqSJ7Z2bzmwyh+FRLbSAaniqSRNydW19dxU0vrmTayN7ccV4BudkqHhI7+naJpIDa+kaun7WEJxcUcfrkPG48Y4IazCXmVEBEktyO3bX88IH5vLtmG1ceO4orjhmlrroSFyogIklszZZKLrr3fYq3V/HHf5vEaZPzwo4kaUQFRCQJuTuzFhRzw5xlZGUaD/3gEA7O7xV2LEkzKiAiSWZbZS0/mbWE55dt4uD8ntxy9iSG9NKsuhJ/KiAiSeSVjzbzX08sYWdVHdedNJYfHDmczAy1d0g4VEBEksC2ylr+37PLeXx+EWMHdOWBi6ey/8CkG4olKUYFRCSBNTY6Mws3cOPzH7Grup4fTR/BlceOUhddSQgqICIJallJOT97eikL1+9gan4vfn3aAYwZ0DXsWCKfUgERSTClO6v5yyureOjddfTMzebmsybyrSl5GtshCUcFRCRBlFfVcdvcT7jnrbXUNTTy3UP245rjx9A9VxMhSmJSAREJ2e7aeu6ft46/vfYJ5VV1zJg4iKuOG01+n85hRxNpkQqISEhWbKrg4XfXMWtBMRU19Rw1pi/XnDCG8YO6hx1NJCoqICJxVF3XwPNLN/HQu+t4f+12sjMzOHnCAM49LJ+D9usZdjyRfaICIhIHa7ZU8sh763m8cAPbd9exX+9cfnLyWM48aAi9tFKgJCkVEJEYqa5r4NWPSnno3fW8uWoLmRnG8eP6851DhjJtRB8yNIJckpwKiEg7WrulktdWlPLayjLeWb2V6rpGBnXP4erjRnP2wUPo3y0n7Igi7UYFROQrqK5rYN7qrcxdUcZrK0pZu3U3APm9c/n2wUM5amw/jhjZR/NVSUpSARHZB42NzorNFcz7ZCtzg6uMmvpGOmZlcNiI3lxweD7Tx/RTF1xJCyogIi2oqW/go40VvLdmG++u2cb7a7dRXlUHwLA+nTln6lCmj+nLocN7k9NB81NJelEBEQlU1TbwSdkulhaXs7i4nCVF5Xy0aSd1DQ5EbkudOH4AU4f1YuqwXlqDQ9KeCoiklZr6Bkp2VFO8vYriHbtZXVbJx6W7WFW6iw3bd+ORWkHXnCwOHNydi48YzoS87hTk91QDuMgXqIBISqmsqad4RxXF26so2lFF0fbdQbGoomh7FWUVNZ87Pjszg+F9O3Pg4O6cMWUwI/t1YfygbgztlatutiKtUAGRhOTu1NQ3sru2gaq6BqpqI4+Kmjq2VdayvbKWrZW1bKuspXRnDUU7IoVi++66z/2cDpnGoB6dyOvRiemj+5LXsxODe+aS16MTg3t2YmD3HLIyM0I6S5HkpgISZ2UVNVw1cxGlO2taPziNNLhTVdtAdV1QMOoaPr2d1JLunTrQt2tH8np0YuLgHuT17PRpccjrkUu/rh11JSESIyogcVRT38C/PzifZSXlfH10Xwz9w7ZHRgbkdNQD90wAAAZaSURBVMikU/DIzc4kJzuT3A6ZdMrOJKdDJrnZWXTOzqRXl2x6dc6mZ242HXT1IBIaFZA4cXd+8fQy5q/bzl+/M5lTDhwUdiQRka9Ev77FyX1vr+Wxwg1cdvRIFQ8RSQkqIHHw5sdb+PU/l3PcuP78x7Gjw44jItIuVEBibO2WSi59eAEj+3bhD/82SQ26IpIykroNxMxOBP4EZAJ3uvuNsfhzrn1iMQvWb2/TZ8t21ZBhcOf5BXTpmNR/3SIin5O0/6KZWSZwK3AcUAS8b2Zz3P3D9v6z8np2oqKmrvUDmzF2YDcumpavaS9EJOUkbQEBpgKr3H01gJk9CpwKtHsBufyYUe39I0VEkl4yt4HkARuabBcF+z7HzC4xs0IzKywrK4tbOBGRVJfMBaS51ugvjV1299vdvcDdC/r27RuHWCIi6SGZC0gRMKTJ9mCgJKQsIiJpJ5kLyPvAKDMbZmbZwLeBOSFnEhFJG0nbiO7u9Wb2Y+AFIt1473b3ZSHHEhFJG0lbQADc/Vng2bBziIiko2S+hSUiIiFSARERkTYxj2bVnhRhZmXAujZ+vA+wpR3jhEHnkBh0DokhFc4B4nMe+7n7l8ZBpFUB+SrMrNDdC8LO8VXoHBKDziExpMI5QLjnoVtYIiLSJiogIiLSJiog0bs97ADtQOeQGHQOiSEVzgFCPA+1gYiISJvoCkRERNpEBURERNpEBSQKZnaima0ws1Vmdl3YeaJhZnebWamZLW2yr5eZvWRmHwfPPcPM2BozG2Jmr5rZcjNbZmZXBPuT5jzMLMfM3jOzD4Jz+FWwf5iZvRucw2PBhKAJzcwyzWyhmf0j2E6qczCztWa2xMwWmVlhsC9pvksAZtbDzJ4ws4+C/y8OC/McVEBa0WTp3JOAccA5ZjYu3FRRuRc48Qv7rgNedvdRwMvBdiKrB6529/2BQ4FLg7/7ZDqPGuBod58ITAJONLNDgd8CfwjOYTtwcYgZo3UFsLzJdjKew1HuPqnJuIlk+i4B/Al43t3HAhOJ/PcI7xzcXY8WHsBhwAtNtq8Hrg87V5TZ84GlTbZXAAOD1wOBFWFn3MfzmQ0cl6znAeQCC4BDiIwczgr2f+47logPIuvtvAwcDfyDyIJuyXYOa4E+X9iXNN8loBuwhqDzUyKcg65AWhfV0rlJor+7bwQInvuFnCdqZpYPTAbeJcnOI7j1swgoBV4CPgF2uHt9cEgyfKf+CPwX0Bhs9yb5zsGBF81svpldEuxLpu/ScKAMuCe4lXinmXUmxHNQAWldVEvnSuyYWRfgSeBKd98Zdp595e4N7j6JyG/xU4H9mzssvqmiZ2anAKXuPr/p7mYOTdhzCExz9ylEbkdfamZfCzvQPsoCpgB/c/fJQCUh33JTAWldKi2du9nMBgIEz6Uh52mVmXUgUjwecvdZwe6kOw8Ad98BvEakPaeHme1ZjyfRv1PTgBlmthZ4lMhtrD+SXOeAu5cEz6XAU0SKeTJ9l4qAInd/N9h+gkhBCe0cVEBal0pL584Bzg9en0+kTSFhmZkBdwHL3f2WJm8lzXmYWV8z6xG87gQcS6Th81XgzOCwhD4Hd7/e3Qe7ez6R7/8r7v5dkugczKyzmXXd8xo4HlhKEn2X3H0TsMHMxgS7jgE+JMRz0Ej0KJjZyUR+49qzdO5vQo7UKjN7BJhOZKrnzcAvgaeBmcBQYD1wlrtvCytja8zsCOANYAmf3Xv/CZF2kKQ4DzM7ELiPyHcnA5jp7v9tZsOJ/DbfC1gIfM/da8JLGh0zmw5c4+6nJNM5BFmfCjazgIfd/Tdm1psk+S4BmNkk4E4gG1gNXEjwvSKEc1ABERGRNtEtLBERaRMVEBERaRMVEBERaRMVEBERaRMVEBERaRMVEJEEYWavmVlB60eKJAYVEBERaRMVEJE4M7P8YC2HO4I1Ql4MRqnveT/DzO4zs/8JM6dIa1RARMIxCrjV3ccDO4Azgv1ZwEPASnf/WVjhRKKhAiISjjXuvih4PZ/I2i0AtxFZwyXhp8sRUQERCUfTOaMaiFx5ALwNHGVmOfGPJLJvVEBEEstdwLPA402mShdJSCogIgkmmLp+AfCAmen/UUlYmo1XRETaRL/diIhIm6iAiIhIm6iAiIhIm6iAiIhIm6iAiIhIm6iAiIhIm6iAiIhIm/x/MxmtDmnq2FEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "old_vgrid.dz.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'depth [m]')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAHgCAYAAAALws5bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhU1ZnH8e/b1QvNJoKA2kBARRKjIkrAZTTuqFFBo4mZSTSbZhI1GpcoYxITE0cNcU2MjkZHTYxLEoJOXAjuSwRBm80FaEGERllkFRro5Z0/6jYWTVV3Qde9t6rr93meeuw6dav6vSn95fS5555j7o6IiISrJO4CRESKgcJWRCQCClsRkQgobEVEIqCwFRGJgMJWRCQCpXEXEJdddtnFBw4cGHcZItLBvPHGGyvcvXfL9qIN24EDBzJt2rS4yxCRDsbMFqZr1zCCiEgEFLYiIhFQ2IqIREBhKyISAYWtiEgEFLYiIhFQ2IqIREBhKyISAYWtiEgEFLYiIhFQ2IqIREBhKyISAYWtiEgEFLYiIhHoMGFrZieY2RwzqzGzK+OuR0QkVYdYz9bMEsDtwHHAYmCqmT3u7m+397OPu+kF5i1bv+X54D5dmHTJke39WBEpMh0ibIERQI27zwcws4eB0UC7wrZl0ALMW7ae4256YUvg/mTCLB6asohGdxJmfG1kf341Zr/2/FoR6YA6SthWAYtSni8GRrb3Q1sGbWr7y/OW8+CUhTw9e+mW9kZ3/jT5AwAFrohspaOM2VqaNt/mILPzzGyamU1bvnx5u37hN+55faugTfVgELgiIs06StguBvqnPO8HLGl5kLvf5e7D3X14797b7Me2XR793iEZX9sm5UWk6HWUsJ0KDDazQWZWDpwFPN7eDx3cp0vG9hGDerb340WkiHSIsHX3BuACYCLwDvCou7/V3s+ddMmR2wSuZiOIyI7oKBfIcPcngSdz/bkKVhHJhQ7RsxURyXcKWxGRCChsRUQioLAVEYmAwlZEJAIKWxGRCChsQ/Ifd78WdwkikkcUtiF59b2V/GTCrLjLEJE8obBth507l7X6+p8mf8CE6tqIqhGRfKawbYerT/l8m8eMHT8zgkpEJN8pbNthzLAqDtuz9QVp6uqbGHbNP9XDFSlyCtt2evDcQ9oM3FUb6vnRI9M1hitSxBS2OfDguYfQpTzR6jFOclFx9XBFipPCNkeuPa3tbXAcuPTRGQpckSKksM2RMcOq+PrBA9Luz5Oq0Z2LH5mucVyRIqOwzaFfjdmPm796AD0qW58SBslx3LHjZylwRYqEwjbHxgyrYvrVx2fVy62rb9SwgkiRMPfi3J5w+PDhPm3atFB/x4TqWi59dAaNbfxvbCTHc6t6VHL5qCGMGVYVal0iEh4ze8Pdh7dsV882RGOGVXHjV4ZSWdb2TAWA2tV1GloQ6aAUtiEbM6yK607fL6txXEgOLVz8yHQOu/45ha5IB6KwjUDzOO4tXz2AhLU1kpukXq5Ix6KwjVC2wwrN6uobGTdxTshViUgUFLYRax5WqOpRCdDmjIXa1XUMuvIJDSuIFLjSuAsoRmOGVW2ZcTChupZxE+dQu7ou4/HOp8MKze8XkcKiqV95YkJ1LWPHz6KuvrHV4xJmNLmzu6aJieSlTFO/1LPNE82hOW7iHJasriPT/wU2z9lVT1eksGjMNo+MGVbFq1cezYLrv7RlTLc1uoAmUjgUtnnq8lFDspq1sGR1HROqazns+ud0IU0kj2kYIU+1HFYoMUt7268Dl/1lBg1NGl4QyWcK2zzWctZCywtoFaUlNLlT37h1CDcPLyhsRfKHhhEKROr8XCO5aM0NX96fhsb0l9I0vCCSXzT1q8Addv1zaefolhiUJkrY3NC0pa2yLMF1p++nHq9IiLTqVweV7kJaWcJwZ6ugBc1eEImTwrbApRteGHfG0IzHL2nlTjURCY8ukHUAqRfSmmW6BbhTWYIFK9YzY9HqLTMddDeaSPgUth3U5aOGbDN7obTEqG9s5OjfvLDVVDJNFxMJn4YROqh0wwu/OXMo/xp7DJ0rEtvM2dV4rki41LPtwNINLwBs2JR+sRuN54qERz3bIrR7hnUXHPjphNl8uEahK5Jr6tkWoXTjuRWlJRw4oAcPvf4Bj0xdxNdG9GdQ7y7c/dICXUQTyQGFbRFque5CapAuWrmB379QwwOvLdxqmUddRBNpH91BJmmN/O9nWLp20zbtVT0qefXKo2OoSKQw6A4y2S7L0gQtJHu4rW3hIyLpKWwlrUwX0QCOHPc8Y8fPZNHKDRFWJFLYNGYraaW7iFZZluDHJwxh/vL1PDJ1EX+ZtpjTD6xi777d+N9X39eFNJFWKGwlrdYuogGcf9Re3Pnie/zxtfdJXeVRF9JE0tMFMmmXg//7WT5au3Gbdl1Ik2KlC2QSiqVpghaSPdyaZZ9EXI1I/lLYSrtkupBmwHE3v8iFD1Uzb+m6aIsSyUMKW2mXdIuXV5Yl+OWYz/OfX9yT595ZyvG3vMT5f36TOR8pdKV46QKZtEtbF9LOPXwP7nllPve9+j5PzPyQk/bblX2rduLByR9o9oIUFV0gk0is3rCZe15ZwF0vvsemFptUam806Uh0gUxi1aNzOZceP4SeXSu2eU1r6UoxUNhKpD5ak3n2wuzaNRFXIxIdha1EqrXZCyf/9hW+e/80ha50SApbiVSm2QvXnr4vPzp2b15f8HEQulOZtVihKx2HZiNIpNqavfCtfxvIfa++zx9ens8pv3uFYz7bh4uOHcz+/XowobpWOwJLwdJsBMlLazfWc/+r7/OHVxawpq6efXbrxnvL17OpoWnLMZrFIPlIsxGkoHTvVMaFxwzmlSuO4rLj9+adj9ZtFbSgWQxSWBS2kte6dSrjgqMHQ4Y/wLQjsBQKha0UhEyzGDqVlWjBGykIClspCOlmMZSWGA1NzvE3v8ilj87gg4+1c4TkL81GkIKQaRbDvw3ehTtfeI8/Tl7IY9NrOXN4fy48eq9Wt/URiYNmI0iHsHTtRm5/voaHXv8Aw/j3kQP4wVF70qdbp7hLkyKTV7MRzOxMM3vLzJrMbHiL18aaWY2ZzTGzUSntJwRtNWZ2ZUr7IDObYmbzzOwRMyuP8lwkP/Tt3olrRu/L85cdyWnDqvjj5IUc8evnue6pd1i1fjMTqms57PrnGHTlExx2/XNMqK6Nu2QpMrH0bM3sc0AT8D/AZe4+LWjfB3gIGAHsDjwD7B28bS5wHLAYmAp8zd3fNrNHgfHu/rCZ3QnMcPc72qpBPduObcGK9dz6zFwem7GE8hKj0aGh6dN/1zVHV8KSVz1bd3/H3dNNkBwNPOzum9x9AVBDMnhHADXuPt/dNwMPA6PNzICjgb8G778fGBP+GUi+G7RLF245axgTLz4CCy6kpdIcXYlavs1GqAIWpTxfHLRlau8FrHb3hhbtIgDs3bcbm+qb0r6mOboSpdDC1syeMbPZaR6jW3tbmjbfgfZMNZ1nZtPMbNry5ctbPwHpMDLNTCgpMSZU19LUVJwXiSVaoYWtux/r7vumeTzWytsWA/1TnvcDlrTSvgLoYWalLdoz1XSXuw939+G9e/fekdOSApRujm55ooS+3Sq4+JHpnHTbyzz7zlKKdWaORCPfhhEeB84yswozGwQMBl4neUFscDDzoBw4C3jck/91PA+cEbz/HKC1MJciNGZYFdedvh9VPSoxoKpHJb8+Y39eueJobj3rAOrqG/nO/dM4487XmDL/47jLlQ4qrtkIpwG/BXoDq4Hp7j4qeO0q4NtAA3Cxuz8VtJ8E3AIkgHvd/dqgfQ+SF8x6AtXA1919U1s1aDaCNKtvbOLRaYu47dl5LF27iSOH9ObyUUP4/O47xV2aFKBMsxF0U4NIoG5zI/e/9j53vPAea+rqOWXo7lxy3N4M2qVL3KVJAVHYtqCwlUzW1NVz90vzueeVBWxubOIrw/tz0TGDmTz/Yy1eLm1S2LagsJW2LFu3kdufq+HPr39AU5NjZroxQtqUVzc1iBSCPt068YvR+/LcpUdSXprQjRHSLgpbkTb079mZjfWNaV/TjRGSLYWtSBYy3RhRXlrCux+tjbgaKUQKW5EspLsxoixhlBicdOvL/GTCLFau3xxTdVIItHi4SBYyLV5+5JDe3PLMPP44eSGPT1/CxcfuzTcO+QxlCfVjZGuajSCSA/OWruOaf7zNy/NWsGfvLvz05H04ckifuMuSGGg2gkiIBvftxgPfHsE95wynyeGb/zuVb/3v67y3XJtRSpLCViRHzIxjPteXiRcfwVUnfY5p769i1M0v8ct/vM2auvq4y5OYaRhBJCQrPtnEjf+cw8NTF7Fz53IuPX5vKksT3Dhpru5C68B0B1kLCluJyuzaNVzzj7d5fcFKjK0XXNZdaB2PxmxFYrJv1U48ct7B9Oxcts3K9roLrXgobEUiYGas2pB+3FZ3oRUHha1IRDLdhVaaMGbXrom4GomawlYkIpnuQqsoLWH07a9y7RNvs2FzQ4Z3S6FT2IpEJN32POPOGMqrVxzDV4b34+6XF3DcTS/x/JxlcZcqIdBsBJE88fqClYwdP5P3lq/nlKG787OT96F3t4q4y5LtpNkIInluxKCePHnR4Vx87GAmzv6IY258gYeDhcul8ClsRfJIRWmCi4/dmycvOpzP7tadK8fP4qy7JlOzbF3cpUk7KWxF8tBefbry8LkHc8OX92PO0nWceOvL3DxpLpsa0i9iLvlPYSuSp0pKjK9+YQDPXPJFTtx3N259dh4n3voyk+d/HHdpsgN0gUykQLwwZxk/mTCbxavq+Orw/gztvxO3P/+e1lnIM1oboQWFrRSiDZsbuPWZedz10vxtbv3VOgv5QbMRRDqAzuWljD3pc+ySZkqY1lnIbwpbkQK0Yt2mtO1aZyF/KWxFClCmdRa6VpRS39gUcTWSDYWtSAFKt85Cwox1mxr4yv+8xqKVG2KqTDJR2IoUoHTrLNz4laHc9rVh1Cz9hJNufZnHptfGXaak0GwEkQ5m0coNXPzIdN5YuIrTD6zimtH70rWiNO6yioZmI4gUif49O/PIeQfzw2MGM6G6lpNve5kZi1bHXVbRU9iKdECliRIuOW5vHjr3YDY3NPHlO/7FnS++p0VtYqSwFenARu7Ri6cuOoLj9unL9U+9yzfuncLStRvjLqsoKWxFOridOpfx+/84kOtP3483F67mhFte4pm3l8ZdVtHRqLlIETAzzhoxgOEDe/LDh6r57gPTOPuQz7Bf1U7c8sw8ra8QAYWtSBHZq09X/n7+ofz66Tnc88oCDLassVC7uo6x42cBKHBDoGEEkSJTUZrgpyfvQ68u5dssZqP1FcKjsBUpUivXb07brvUVwqGwFSlSmdZX2HWnThFXUhwUtiJFKt36CgBNTa61FUKgsBUpUunWV/jBUXtSV9/ImNtf5Y2Fq+IusUPR2ggispX3ln/Ct++byodrNnLjmUM5ZejucZdUULQ2gohkZc/eXfn7Dw5jaL+duPChan777DyKtVOWSwpbEdlGzy7l/Om7IzltWBU3TprLpX+ZoW3U20k3NYhIWhWlCW76ylAG9urCzc/MZfGqOv7n6wexc5fyuEsrSOrZikhGZsZFxw7m1rMOYPoHqznt968yf/kncZdVkBS2ItKm0QdU8edzR7J2YwOn3/EvJs//OO6SCo7CVkSyMnxgTyb84DB6dSnnG/dM4W9vLI67pIKiMVsRydqAXp0Z/4PD+MGDb3DpX2awYMV69tylC7+ZNFcrh7VBYSsi22WnyjLu+9YIfjphNr97voaEQWMwM0wrh2WmYQQR2W5liRKuO30/uncq3RK0zbRyWHoKWxHZIWbGuo0NaV/TymHbUtiKyA7LtHJYpvZiprAVkR2WbuWwEoMfHTs4poryl8JWRHZYy5XDenQuo8nhuTnLaNS26VvRbAQRaZcxw6q2mnnwh5fn86sn3qFrxUxu+PL+mFmM1eUPha2I5NR3D9+DtXX13PZcDd07lXHVlz6nwEVhKyIh+NFxe7N2YwN/eGUBO1WWceExGsNV2IpIzpkZPzt5H9ZurOfGSXPpXlnGOYcOjLusWClsRSQUJSXGr7+8P+s2NnD142/RvbKU04b1i7us2Gg2goiEpjRRwm+/NoxD9+zFZX+ZyT/f+ijukmKjsBWRUHUqS3DX2cPZt2onLniomn+9tyLukmKhsBWR0HWtKOX+b32Bgb06c+7905i+aHXcJUVOu+uKSGSWrt3ImXe+xtqN9XzviD340+QPOtzSjNpdV0Ri17d7J/70nZE0NjZxw9NzqF1dh/Pp0owTqmvjLjE0sYStmY0zs3fNbKaZ/d3MeqS8NtbMasxsjpmNSmk/IWirMbMrU9oHmdkUM5tnZo+YmXajE8ljA3p1prJ824lQHX1pxrh6tpOAfd19f2AuMBbAzPYBzgI+D5wA/N7MEmaWAG4HTgT2Ab4WHAtwA3Czuw8GVgHfifRMRGS7LV+3KW17R16aMZawdfd/unvzQpiTgebJd6OBh919k7svAGqAEcGjxt3nu/tm4GFgtCXvATwa+Gvw/vuBMVGdh4jsmGJcmjEfxmy/DTwV/FwFLEp5bXHQlqm9F7A6Jbib20Ukj6VbmrE8UcLlo4bEVFH4QruDzMyeAXZN89JV7v5YcMxVQAPwYPPb0hzvpP8/BW/l+Ew1nQecBzBgwICMtYtIuJpnHYybOIclq+tIlBjlpcahe/WKubLwhBa27n5sa6+b2TnAycAx/un8s8VA/5TD+gFLgp/Tta8AephZadC7TT0+XU13AXdBcupX9mcjIrmWujTj3KXrOPV3r3DJIzN44NsjKCnpeKuExTUb4QTgCuBUd9+Q8tLjwFlmVmFmg4DBwOvAVGBwMPOgnORFtMeDkH4eOCN4/znAY1Gdh4jkxt59u/HzUz7PKzUruOPF9+IuJxRxjdn+DugGTDKz6WZ2J4C7vwU8CrwNPA2c7+6NQa/1AmAi8A7waHAsJEP7EjOrITmGe0+0pyIiufDVL/Tn5P1346ZJc3lj4cq4y8k53UEmInlj7cZ6Tr7tFRqbnCd++G/06Fx40+Z1B5mI5L3uncr47deGsWzdRn7815l0pM6gwlZE8srQ/j244oTP8s+3l/LHyQvjLidnFLYikne+fdggjhrSm1/94x3eWrIm7nJyQmErInmnpMT4zZlD2blLGRf+uZr1mxraflOeU9iKSF7q1bWCW746jPc/Xs9PH5sddzntpj3IRCRvHbJnLy48ejC3PjuPTqUlvDh3RcGuf6uerYjktQuP3os9dunCn19fVNDr3ypsRSSvlSZK2LB52zHbQlv/VmErInlv6drCX/9WYSsiea8jrH+rsBWRvJdu/dvKskRBrX+r2QgikvdS17+tXV2HAf/1pc8W1GwEha2IFITm9W/fW/4Jx930IguWb2j7TXlEwwgiUlD27N2VMw7qx5+mLNQFMhGRMP3wmMG4O799bl7cpWRNYSsiBaffzp359xEDeHTaYt5fsT7ucrKisBWRgnT+0XtRljBueWZu3KVkRWErIgWpT7dOfPPQQTw2YwlzPloXdzltUtiKSMH6zy/uQdfyUm78Z/7ftquwFZGC1aNzOd89fA/++fZSZixaHXc5rVLYikhB+87hg+jZpZzf5HnvVmErIgWta0Up3//inrw8bwWT538cdzkZZbyDzMxmZvH+5e5+TA7rERHZbt845DP89rm5nH3P69Q3NuXl4uKt3a6bAE5q5XUDHs9tOSIi2+/p2R9RV99EfWNy6/PmxcWBvAnc1sL2e+7e6j7CZvaDHNcjIrLdxk2csyVomzUvLp4vYZtxzNbdX2nrzdkcIyIStkxrJOTT2gltXiAzs5PNrNrMVprZWjNbZ2ZroyhORCQbhbC4eDazEW4BzgF6uXt3d+/m7t1DrktEJGvpFhfvVFaSV4uLZ7Oe7SJgtrt7m0eKiMQgdXHxJcEOvKcO3T1vxmshu7D9MfCkmb0IbNl1zd1vCq0qEZHt1Ly4uLsz6paXmLv0k7hL2ko2wwjXAhuATkC3lIeISN4xM74yvD/TF61m7tL8WaAmm55tT3c/PvRKRERy5LRhVdzw9Ls8OnURPzl5n7jLAbLr2T5jZgpbESkYvbpWcOzn+jK+upbNDU1xlwNkF7bnA0+bWZ2mfolIofjKF/qzcv1mnn1nadylAFmEbTDVq8TdKzX1S0QKxRGDe7Nr9048Om1R3KUArYStme3a1puzOUZEJA6JEuOMg/rx4tzlfLRmY9zltNqzfTKL92dzjIhILM4c3o8mh7+9uTjuUloN26HBGG2mxzqgb1SFiohsr8/06sLBe/Tk0WmLaGqK976s1haiSQRjtJke3dw9f27PEBFJ46tf6M/CjzcwZcHKWOvQTg0i0qGduO9uVJQa375vKoOufILDrn+OCdW1kdeRzU0NIiIF6+nZH9HQBI1NjUB8C4urZysiHdq4iXNobEq/sHiUsurZmlmC5MWwLce7+wdhFSUikiv5srB4m2FrZhcCVwNLgeb73hzYP8S6RERyYvceldSmCdaoFxbPZhjhImCIu3/e3fcLHgpaESkI6RYWryxLRL6weLaLh68JuxARkTA0XwS74el3+XDNRrp1KuWXo/eNfGHxjGFrZpcEP84HXjCzJ9Di4SJSgJoXFv/SbS/TvVNZLDs4tDaM0LxI+AfAJKA8pa1r+KWJiOTWyEG9ePODVWxqaIz8d2fs2br7LwDM7Ex3/0vqa2Z2ZtiFiYjk2sg9enLvqwuYuXgNXxjYM9Lfnc0FsrFZtomI5LURQcBOmf9x5L+7tTHbE4GTgCozuy3lpe5AQ9iFiYjk2s5dyvnsrt2YsmAlF0T8u1vr2S4BpgEbgTdSHo8Do8IvTUQk90YO6skbC1dR3xjtdjmtrfo1w93vB/YCHgKqgTeBf7j7qojqExHJqZF79GLD5kZm1UY7ozWbMdvjgPeA24DfATXBEIOISMEZMah53DbaJRezCdubgKPc/Uh3/yJwFHBzuGWJiIRjl64V7NWnK1MWRHuRLJuwXebuNSnP5wPLQqpHRCR0Iwf1ZNr7q2iIcNw2m7B9y8yeNLNvmtk5wP8BU83sdDM7PeT6RERybuQevfhkUwNvf7g2st+ZzdoInUiu+PXF4PlyoCdwCsnVv8aHU5qISDhWr0+uPHDq716lqkcll48aEvotvG2Grbt/K9QKREQiNKG6luue+nTh8Kh2bmhzGMHM9jazZ81sdvB8fzP7SWgViYiEaNzEOdTVb702QhQ7N2QzZns3ydtz6wHcfSZwVphFiYiEJa6dG7IJ287u/nqLNt2uKyIFKdMODWHv3JBN2K4wsz1JXgzDzM4APgy1KhGRkMS1c0M2sxHOB+4CPmtmtcAC4OuhViUiEpLmi2DjJs6hdnUdlWUlXHf6fnkxG2E+cKyZdQFK3H1dqBWJiISseeeGcx+YxoIV6yPZuSGbbXFatgPt2xbHzH4JjCa5W+8y4JvuvsSSH34ryaUdNwTtbwbvOQdongXxq2CRHMzsIOA+oBJ4ErjI3bfeJF5EJI0hfbvx/LvL2NTQSEVpou03tEM22+IMB74PVAWP/wT2aefvHefu+7v7AcA/gJ8F7ScCg4PHecAdAGbWk+R26iOBEcDVZrZz8J47gmOb33dCO2sTkSKx967daGhyFqxYH/rvam2JxV8EW+PsAhzo7pe6+6XAQUC/9vxSd0+9R64LwcU3kr3dBzxpMtDDzHYjuX7uJHdfGSzvOAk4IXitu7u/FvRmHwDGtKc2ESkeQ/p2A2DOR+GPjmZzgWwAsDnl+WZgYHt/sZldC5xNcpv0o4LmKpJbpzdbzKc96kzti9O0Z/qd55HsBTNgwID2nYCIFLxBu3ShtMSYuzT8sM1m6tcfgdfN7OdmdjUwBbi/rTeZ2TNmNjvNYzSAu1/l7v2BB2HLDhWW5qN8B9rTcve73H24uw/v3bt3W6cgIh1ceWkJe/TuwpyPPgn9d2UzG+FaM3sKODxo+pa7V2fxvmOzrOHPwBMkx2QXA/1TXutHcnuexcCRLdpfCNr7pTleRCQrg/t2Y9bi8HdtyKZni7u/6e63Bo82g7YtZjY45empwLvBz48DZ1vSwcAad/8QmAgcb2Y7BxfGjgcmBq+tM7ODg5kMZwOPtbc+ESkeQ/p2Y9GqDWzYHO6NsdmM2YbhejMbQnLq10KSMxwgOXXrJKCG5NSvbwG4+8pgutjU4Lhr3L15T4vv8+nUr6eCh4hIVvbu2w13qFn2Cfv36xHa74klbN39yxnaneQda+leuxe4N037NGDfnBYoIkVjyK6fzkgIM2yzGkYQEemoBvTsTEVpSegzEhS2IlLU/m/GEprcufvlBRx2/XNMqK4N5fcobEWkaE2ormXs+FnUNyZnjDbv2hBG4CpsRaRoRblrg8JWRIpWlLs2KGxFpGhFuWuDwlZEilaUuzbEdVODiEjsmhcN//XEd1myeiPdOpXyy9H7hrKYuHq2IlLUxgyr4l9XHkOvLuWcMnT30HZtUNiKiAC9u1WwbO2m0D5fYSsiAvTp3oll6zaG9vkKWxERoK96tiIi4evTvYLln2yiqSmc/WIVtiIiQJ9unWhscj5ev7ntg3eAwlZEBOjbvQIgtHFbha2ICNC7WycAlq0LZ9xWYSsiAvTpFvRs16pnKyISmj7NwwghzUhQ2IqIABWlCXp0LmOpxmxFRMIzobqWTzY28KfJH4SyY4PCVkSKXvOODQ1N4e3YoLAVkaIXxY4NClsRKXpR7NigsBWRohfFjg0KWxEpelHs2KCdGkSk6DUvGH7N/73Fyg319OlWwX+d9LmcLiSusBURIRm4fbpX8O93T+HWs4ZxyJ69cvr5GkYQEQn0qCwHYE1d7lf+UtiKiAR6dC4DYPWG+px/tsJWRCSwJWzrFLYiIqGpLEtQnihRz1ZEJExmxk6dyzRmKyISth6VZerZioiErUdnha2ISOh2qizXBTIRkbD16FzGmg0asxURCVWPyjL1bEVEwjShupZHpi1iw+ZGDr3u2ZwuHq61EURE+HS3huZFxJes2cjY8bMAcrIgjXq2IiKEv1uDwlZEhPB3a1DYiogQ/m4NCqrUfxwAAA1BSURBVFsREcLfrUEXyERE+PQi2LVPvsPydZvo2aWcn528T852a1DPVkQkMGZYFQ+dOxKAn5/6+Zxui6OwFRFJUZ5IDiVsbmjK6ecqbEVEUpSXJmNRYSsiEqJPw7axjSO3j8JWRCRFc9huUs9WRCQ8FRpGEBEJX2mJYQabGxW2IiKhMTPKEyXq2YqIhK28tERjtiIiYasoLdEwgohI2MoTJWyqV9iKiISqPISerRaiERFJMaG6lsWr6nj/4w28uXAVl48aop0aRERyqXlrnIYmB6B2dR1jx8/KyV5kClsRkUCYW+MobEVEAmFujaOwFREJhLk1jsJWRCQQ5tY4mo0gIhJonnUwdvxM6uqbqOpRmbPZCApbEZEUY4ZV8dK85UyZv5JXrzw6Z58b6zCCmV1mZm5muwTPzcxuM7MaM5tpZgemHHuOmc0LHuektB9kZrOC99xmZhbHuYhIx9Ghbtc1s/7AccAHKc0nAoODx3nAHcGxPYGrgZHACOBqM9s5eM8dwbHN7zshivpFpOPqaKt+3Qz8GPCUttHAA540GehhZrsBo4BJ7r7S3VcBk4ATgte6u/tr7u7AA8CYaE9DRDqairJExwhbMzsVqHX3GS1eqgIWpTxfHLS11r44TbuIyA4rTxTQ2ghm9gywa5qXrgL+Czg+3dvStPkOtGeq6TySQw4MGDAg02EiUuTKS0tobHIaGpsoTeSmTxpa2Lr7senazWw/YBAwI7iW1Q9408xGkOyZ9k85vB+wJGg/skX7C0F7vzTHZ6rpLuAugOHDh2cMZREpblt22M1h2EY+jODus9y9j7sPdPeBJAPzQHf/CHgcODuYlXAwsMbdPwQmAseb2c7BhbHjgYnBa+vM7OBgFsLZwGNRn5OIdCzlidxv+phv82yfBE4CaoANwLcA3H2lmf0SmBocd427rwx+/j5wH1AJPBU8RER2WEVZBwzboHfb/LMD52c47l7g3jTt04B9w6pPRIpPc882l/uQaW0EEZEWUsdsc0VhKyLSwvRFqwA49sYXOez657R4uIhIrk2oruXByclp/U7udmtQ2IqIpBg3cc42wwe52K1BYSsikiKs3RoUtiIiKcLarUFhKyKS4vJRQ6go3Toac7Fbg8JWRCTFmGFVXHDUXlueV/Wo5LrT92v3bg2x39QgIpJvjvlcX26cNJc7v34QJ+ybbj2t7aeerYhIC4mS5IKCTZ679aoUtiIiLTSHbUOTwlZEJDRberYKWxGR8JSqZysiEr6SIGwbm7QQjYhIaEq3hG3uPlNhKyLSQkI9WxGR8CWsOWw1ZisiEppEQhfIRERC19yz1U0NIiIhenLWhwD895PvaqcGEZEwTKiu5aePzd7yXDs1iIiEYNzEOWys104NIiKh0k4NIiIR0E4NIiIRuHzUECrLElu15WKnBi0eLiKSonlHhh89Mh0nuVPD5aOGtHunBvVsRURaGDOsim6dSvnmoQN59cqj2x20oLAVEUkrUWK6XVdEJGyJkhIadQeZiEi4EiXQ2KiwFREJVWlJiRaiEREJW6LEtBCNiEjYEiWmnq2ISNgSJabddUVEwpYwo0Hb4oiIhEvzbEVEIlCaUNiKiISuxHSBTEQkVBOqa3l7yVpenrdC2+KIiIRhQnUtY8fPYnNj8uKYtsUREQnBuIlzqKtv3KpN2+KIiOSYtsUREYmAtsUREYmAtsUREYlA864MV/xtJpsamnK2LY7CVkSkhTHDqvjrG4upq2/kb98/NCefqWEEEZEIKGxFRCKgsBURiYDCVkQkAgpbEZEIKGxFRCKgsBURiYDCVkQkAgpbEZEIKGxFRCKgsBURaWFCdS1T31/JGwtXaacGEZEwNO/UsKlBOzWIiIRGOzWIiERAOzWIiERAOzWIiERAOzWIiERAOzWIiEREOzWIiBSoWMLWzH5uZrVmNj14nJTy2lgzqzGzOWY2KqX9hKCtxsyuTGkfZGZTzGyemT1iZuVRn4+ISFvi7Nne7O4HBI8nAcxsH+As4PPACcDvzSxhZgngduBEYB/ga8GxADcEnzUYWAV8J+oTERFpS74NI4wGHnb3Te6+AKgBRgSPGnef7+6bgYeB0WZmwNHAX4P33w+MiaFuEZFWxRm2F5jZTDO718x2DtqqgEUpxywO2jK19wJWu3tDi/a0zOw8M5tmZtOWL1+eq/MQEWlTaGFrZs+Y2ew0j9HAHcCewAHAh8CNzW9L81G+A+1puftd7j7c3Yf37t17u85HRKQ9Qpv65e7HZnOcmd0N/CN4uhjon/JyP2BJ8HO69hVADzMrDXq3qceLiOSNuGYj7Jby9DRgdvDz48BZZlZhZoOAwcDrwFRgcDDzoJzkRbTH3d2B54EzgvefAzwWxTmIiGyPuG5q+LWZHUDyT/73ge8BuPtbZvYo8DbQAJzv7o0AZnYBMBFIAPe6+1vBZ10BPGxmvwKqgXuiPBERkWzEErbu/o1WXrsWuDZN+5PAk2na55OcrSAikrfybeqXiEiHpLAVEWlB2+KIiIRM2+KIiERA2+KIiERA2+KIiERA2+KIiERA2+KIiERA2+KIiERE2+KIiBQoha2ISAQUtiIiEVDYiohEQGErIhIBha2ISAQUtiIiEVDYiohEQGErIhIBha2ISAQUtiIiLWinBhGRkGmnBhGRCGinBhGRCGinBhGRCGinBhGRCGinBhGRCDTvyDBu4hyWrK5jd+3UICISjjHDqtodri1pGEFEJAIKWxGRCChsRUQioLAVEYmAwlZEJAIKWxGRCChsRUQioLAVEYmAwlZEJAIKWxGRCChsRUQioLAVEYmAwlZEJAIKWxGRCChsRUQiYO4edw2xMLPlwMIsDt0FWBFyOVHRueSvjnQ+xX4un3H33i0bizZss2Vm09x9eNx15ILOJX91pPPRuaSnYQQRkQgobEVEIqCwbdtdcReQQzqX/NWRzkfnkobGbEVEIqCerYhIBBS2GZjZCWY2x8xqzOzKuOvZEWb2vpnNMrPpZjYtaOtpZpPMbF7wz53jrjMdM7vXzJaZ2eyUtrS1W9JtwXc108wOjK/ybWU4l5+bWW3w3Uw3s5NSXhsbnMscMxsVT9XpmVl/M3vezN4xs7fM7KKgveC+m1bOJZzvxt31aPEAEsB7wB5AOTAD2CfuunbgPN4HdmnR9mvgyuDnK4Eb4q4zQ+1HAAcCs9uqHTgJeAow4GBgStz1Z3EuPwcuS3PsPsG/bxXAoODfw0Tc55BS327AgcHP3YC5Qc0F9920ci6hfDfq2aY3Aqhx9/nuvhl4GBgdc025Mhq4P/j5fmBMjLVk5O4vAStbNGeqfTTwgCdNBnqY2W7RVNq2DOeSyWjgYXff5O4LgBqS/z7mBXf/0N3fDH5eB7wDVFGA300r55JJu74bhW16VcCilOeLaf1LyFcO/NPM3jCz84K2vu7+IST/ZQP6xFbd9stUe6F+XxcEf1rfmzKcUzDnYmYDgWHAFAr8u2lxLhDCd6OwTc/StBXitI3D3P1A4ETgfDM7Iu6CQlKI39cdwJ7AAcCHwI1Be0Gci5l1Bf4GXOzua1s7NE1bXp1PmnMJ5btR2Ka3GOif8rwfsCSmWnaYuy8J/rkM+DvJP3mWNv8ZF/xzWXwVbrdMtRfc9+XuS9290d2bgLv59M/RvD8XMysjGU4Puvv4oLkgv5t05xLWd6OwTW8qMNjMBplZOXAW8HjMNW0XM+tiZt2afwaOB2aTPI9zgsPOAR6Lp8Idkqn2x4GzgyvfBwNrmv+kzVctxi1PI/ndQPJczjKzCjMbBAwGXo+6vkzMzIB7gHfc/aaUlwruu8l0LqF9N3FfEczXB8mrqHNJXnG8Ku56dqD+PUheOZ0BvNV8DkAv4FlgXvDPnnHXmqH+h0j+CVdPskfxnUy1k/zz7vbgu5oFDI+7/izO5Y9BrTOD/4h3Szn+quBc5gAnxl1/i3P5N5J/Os8EpgePkwrxu2nlXEL5bnQHmYhIBDSMICISAYWtiEgEFLYiIhFQ2IqIREBhKyISAYWtiEgEFLZSdIIl9C7bjmNrzeya7fwdD5rZSjM7Y8eqlI5GYSvStpvd/Wfb8wZ3/w8K7K5DCVdp3AWIRMHMrgLOJrlq03LgDTPbHXgy5bD9gD3cfWErn/NzkmuZ7gbsDVxCcp3WE4Fa4BR3rw/jHKSwqWcrHZ6ZHURyfYthwOnAFyC5UI+7H+DuB5BccORvrQVtij2BL5Fc3/RPwPPuvh9QF7SLbEM9WykGhwN/d/cNAGa21Z/3ZnYY8N3guGw85e71ZjaL5K4eTwfts4CBOalYOhyFrRSLtIuABCs83QOc6u6fZPlZmwDcvcnM6v3TBUaa0H9TkoGGEaQYvAScZmaVwbKTp8CWtUwfBa5w97lxFigdn8JWOjxP7jP1CMkl9P4GvBy8dCjJ8dtfpOykuntMZUoHpyUWRVoRzD74xN1/swPvvQ/4h7v/Ndd1SeFRz1akdZ8A5+3ITQ3AF4GNoVQlBUc9WxGRCKhnKyISAYWtiEgEFLYiIhFQ2IqIREBhKyISgf8HR2ckqRC6IJ0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# interfaces\n",
    "e = np.zeros(len(old_vgrid.dz)+1)\n",
    "e[1::] = old_vgrid.dz.cumsum()\n",
    "eh = 0.5 * (e[0:-1]+e[1::])\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,8))\n",
    "ax.plot(old_vgrid.dz, -eh, 'o-')\n",
    "ax.set_xlabel('dz [m]')\n",
    "ax.set_ylabel('depth [m]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.DataArray &#x27;dz&#x27; (nk: 63)&gt;\n",
       "array([2.50000e+00, 5.00000e+00, 7.50000e+00, 1.00000e+01, 2.00000e+01,\n",
       "       3.00000e+01, 4.00000e+01, 5.00000e+01, 6.00000e+01, 7.00000e+01,\n",
       "       8.00000e+01, 9.00000e+01, 1.00000e+02, 1.10000e+02, 1.20000e+02,\n",
       "       1.30000e+02, 1.40000e+02, 1.50000e+02, 1.60000e+02, 1.70250e+02,\n",
       "       1.80750e+02, 1.91750e+02, 2.03500e+02, 2.16000e+02, 2.29250e+02,\n",
       "       2.43500e+02, 2.59000e+02, 2.76000e+02, 2.94750e+02, 3.15250e+02,\n",
       "       3.38250e+02, 3.63750e+02, 3.92500e+02, 4.25000e+02, 4.61750e+02,\n",
       "       5.03750e+02, 5.51750e+02, 6.07000e+02, 6.70750e+02, 7.44500e+02,\n",
       "       8.29750e+02, 9.28250e+02, 1.04150e+03, 1.17125e+03, 1.31825e+03,\n",
       "       1.48300e+03, 1.66500e+03, 1.86300e+03, 2.07475e+03, 2.29800e+03,\n",
       "       2.52975e+03, 2.76800e+03, 3.01050e+03, 3.25600e+03, 3.50325e+03,\n",
       "       3.75175e+03, 4.00075e+03, 4.25025e+03, 4.50000e+03, 4.75000e+03,\n",
       "       5.00000e+03, 5.25000e+03, 5.50000e+03], dtype=float32)\n",
       "Dimensions without coordinates: nk</pre>"
      ],
      "text/plain": [
       "<xarray.DataArray 'dz' (nk: 63)>\n",
       "array([2.50000e+00, 5.00000e+00, 7.50000e+00, 1.00000e+01, 2.00000e+01,\n",
       "       3.00000e+01, 4.00000e+01, 5.00000e+01, 6.00000e+01, 7.00000e+01,\n",
       "       8.00000e+01, 9.00000e+01, 1.00000e+02, 1.10000e+02, 1.20000e+02,\n",
       "       1.30000e+02, 1.40000e+02, 1.50000e+02, 1.60000e+02, 1.70250e+02,\n",
       "       1.80750e+02, 1.91750e+02, 2.03500e+02, 2.16000e+02, 2.29250e+02,\n",
       "       2.43500e+02, 2.59000e+02, 2.76000e+02, 2.94750e+02, 3.15250e+02,\n",
       "       3.38250e+02, 3.63750e+02, 3.92500e+02, 4.25000e+02, 4.61750e+02,\n",
       "       5.03750e+02, 5.51750e+02, 6.07000e+02, 6.70750e+02, 7.44500e+02,\n",
       "       8.29750e+02, 9.28250e+02, 1.04150e+03, 1.17125e+03, 1.31825e+03,\n",
       "       1.48300e+03, 1.66500e+03, 1.86300e+03, 2.07475e+03, 2.29800e+03,\n",
       "       2.52975e+03, 2.76800e+03, 3.01050e+03, 3.25600e+03, 3.50325e+03,\n",
       "       3.75175e+03, 4.00075e+03, 4.25025e+03, 4.50000e+03, 4.75000e+03,\n",
       "       5.00000e+03, 5.25000e+03, 5.50000e+03], dtype=float32)\n",
       "Dimensions without coordinates: nk"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "old_vgrid.dz.cumsum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total depth is 5500.0 m\n"
     ]
    }
   ],
   "source": [
    "print('Total depth is {} m'.format(old_vgrid.dz.cumsum()[-1].values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This grid does not meet the criteria stated before since 5500 < 6000. MOM6 will inflate the last layer such that the total water column is 6000 m. We want to avoid that when creating the new grid by making sure dz.cumsum() = 6000 m.\n",
    "\n",
    "Below is an example showing that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/glade/u/home/gmarques/miniconda3/envs/analysis/lib/python3.7/site-packages/xarray/coding/times.py:426: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range\n",
      "  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)\n"
     ]
    }
   ],
   "source": [
    "path42 = '/glade/scratch/gmarques/g.c2b6.GJRA.TL319_t061.long_JRA_mct.042/run/'\n",
    "ds_hm = xr.open_dataset(path42+'g.c2b6.GJRA.TL319_t061.long_JRA_mct.042.mom6.hm_0001.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  2.5 ,   2.5 ,   2.5 ,   2.5 ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.  ,  10.  ,  10.25,  10.5 ,  11.  ,  11.75,  12.5 ,\n",
       "        13.25,  14.25,  15.5 ,  17.  ,  18.75,  20.5 ,  23.  ,  25.5 ,\n",
       "        28.75,  32.5 ,  36.75,  42.  ,  48.  ,  55.25,  63.75,  73.75,\n",
       "        85.25,  98.5 , 113.25, 129.75, 147.  , 164.75, 182.  , 198.  ,\n",
       "       211.75, 223.25, 231.75, 238.25, 242.5 , 245.5 , 247.25, 248.5 ,\n",
       "       249.  , 249.5 , 249.75, 250.  , 250.  , 250.  , 750.  ])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_hm = np.diff(ds_hm.zi)\n",
    "dz_hm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shown the last layer is 750 m thick and not 250 m as before. Let's make sure the total depth is indeed 6000m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6000.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_hm.cumsum()[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to modifify the old grid to make sure old_vgrid.dz.cumsum()[-1] = 6000 m. An easy way of meeting this criteria is two add two additional levels at the bottom with dz = 250 m each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "dz_new = np.ones(len(old_vgrid.dz) + 2) * 250.\n",
    "dz_new[0:len(old_vgrid.dz)] = old_vgrid.dz.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.50000e+00, 5.00000e+00, 7.50000e+00, 1.00000e+01, 2.00000e+01,\n",
       "       3.00000e+01, 4.00000e+01, 5.00000e+01, 6.00000e+01, 7.00000e+01,\n",
       "       8.00000e+01, 9.00000e+01, 1.00000e+02, 1.10000e+02, 1.20000e+02,\n",
       "       1.30000e+02, 1.40000e+02, 1.50000e+02, 1.60000e+02, 1.70250e+02,\n",
       "       1.80750e+02, 1.91750e+02, 2.03500e+02, 2.16000e+02, 2.29250e+02,\n",
       "       2.43500e+02, 2.59000e+02, 2.76000e+02, 2.94750e+02, 3.15250e+02,\n",
       "       3.38250e+02, 3.63750e+02, 3.92500e+02, 4.25000e+02, 4.61750e+02,\n",
       "       5.03750e+02, 5.51750e+02, 6.07000e+02, 6.70750e+02, 7.44500e+02,\n",
       "       8.29750e+02, 9.28250e+02, 1.04150e+03, 1.17125e+03, 1.31825e+03,\n",
       "       1.48300e+03, 1.66500e+03, 1.86300e+03, 2.07475e+03, 2.29800e+03,\n",
       "       2.52975e+03, 2.76800e+03, 3.01050e+03, 3.25600e+03, 3.50325e+03,\n",
       "       3.75175e+03, 4.00075e+03, 4.25025e+03, 4.50000e+03, 4.75000e+03,\n",
       "       5.00000e+03, 5.25000e+03, 5.50000e+03, 5.75000e+03, 6.00000e+03])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_new.cumsum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  2.5 ,   2.5 ,   2.5 ,   2.5 ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,\n",
       "        10.  ,  10.  ,  10.  ,  10.25,  10.5 ,  11.  ,  11.75,  12.5 ,\n",
       "        13.25,  14.25,  15.5 ,  17.  ,  18.75,  20.5 ,  23.  ,  25.5 ,\n",
       "        28.75,  32.5 ,  36.75,  42.  ,  48.  ,  55.25,  63.75,  73.75,\n",
       "        85.25,  98.5 , 113.25, 129.75, 147.  , 164.75, 182.  , 198.  ,\n",
       "       211.75, 223.25, 231.75, 238.25, 242.5 , 245.5 , 247.25, 248.5 ,\n",
       "       249.  , 249.5 , 249.75, 250.  , 250.  , 250.  , 250.  , 250.  ,\n",
       "       250.  ])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# new interfaces\n",
    "e = np.zeros(len(dz_new)+1)\n",
    "e[1::] = dz_new.cumsum()\n",
    "eh_new = (0.5 * (e[0:-1]+e[1::]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New grid as 65 layers; old grid has 63 layers\n"
     ]
    }
   ],
   "source": [
    "print('New grid as {} layers; old grid has {} layers'.format(len(dz_new), len(old_vgrid.dz)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Can we represent this profile using an analytical function?\n",
    "\n",
    "The answer is yes! The hyperbolic tangent function does a pretty good job at representing this profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hyperbolic tangent function\n",
    "def func_tanh(z,a,b,c,d):\n",
    "    '''\n",
    "    a = dz transitional\n",
    "    b = dz range\n",
    "    c = maximum slope\n",
    "    d = depth where dz transitional should be applied\n",
    "    '''\n",
    "    return a + b * (np.tanh(c/b * (d-z)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial guesses\n",
    "h_max = 6000. # mximum depth\n",
    "a = 100.    # dz transitional\n",
    "b  = 150.  #(dz_max - dz_min)\n",
    "c = -b/1500.  # maximum slope\n",
    "d = 1000.      # depth where dz_t should be applied\n",
    "z = np.linspace(0,h_max, 63)\n",
    "\n",
    "dz = func_tanh(z,a,b,c,d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.Dataset(data_vars={ 'dz' : (('z'), dz)}, coords={'z': z})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'depth [m]')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAHgCAYAAAALws5bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUVeLG8e+ZSQcSekuABKS3CKFXRQGxoFhQUUFFWLtrL7u6q+vq7rp2RbEBLir8rFhBLBRFpAVBkCK9dwgJ6ef3x51gCEkYyNxJez/PM8+QO/eec0Z8Xk7OPfccY61FRETc5SntBoiIVAYKWxGRIFDYiogEgcJWRCQIFLYiIkGgsBURCYKQ0m5Aaaldu7aNj48v7WaISAWzaNGiPdbaOgWPV9qwjY+PZ+HChaXdDBGpYIwxGws7rmEEEZEgUNiKiASBwlZEJAgq7ZitlI6srCy2bNlCenp6aTdFpEQiIiKIi4sjNDTUr/MVthJUW7ZsoVq1asTHx2OMKe3miJwSay179+5ly5YtJCQk+HWNhhEkqNLT06lVq5aCVso1Ywy1atU6qd/QFLYSdApaqQhO9v9jha1IEYYMGcKBAweKPefhhx9m5syZp1T+999/z3nnnXdK10r5ozFbkQKstVhr+eKLL0547qOPPhqEFklFoJ6tVEpPP/007dq1o127djz77LNs2LCB1q1bc9NNN9GpUyc2b95MfHw8e/bsAeCxxx6jVatWnH322VxxxRU89dRTAIwaNYr3338fcJ5KfOSRR+jUqRPt27fnt99+A+Dnn3+mZ8+enH766fTs2ZNVq1aVzpeWUqWerZSav3/6Kyu2HQpomW0aRvPI+W2LPWfRokW89dZbzJ8/H2st3bp1o1+/fqxatYq33nqLl19++ZjzFy5cyAcffMCSJUvIzs6mU6dOdO7cudCya9euzeLFi3n55Zd56qmneP3112nVqhWzZ88mJCSEmTNn8uCDD/LBBx8E7DtL+aCwlUpn7ty5XHTRRVSpUgWAYcOGMWfOHJo0aUL37t0LPX/o0KFERkYCcP755xdZ9rBhwwDo3LkzH374IQAHDx5k5MiRrFmzBmMMWVlZgf5KUg4obKXUnKgH6paiNjnNC19/zy9MeHg4AF6vl+zsbAD++te/csYZZ/DRRx+xYcMG+vfvf3INlgpBY7ZS6fTt25ePP/6YtLQ0UlNT+eijj+jTp0+R5/fu3ZtPP/2U9PR0Dh8+zOeff35S9R08eJDY2FgAJkyYUJKmSzlWYcLWGDPYGLPKGLPWGHN/abdHyq5OnToxatQounbtSrdu3Rg9ejQ1atQo8vwuXbpwwQUX0LFjR4YNG0ZSUhIxMTF+13fvvffywAMP0KtXL3JycgLxFaQcMifzK1JZZYzxAquBs4EtwALgCmvtiqKuSUpKslrPNvhWrlxJ69atS7sZJ+3w4cNUrVqVtLQ0+vbty/jx4+nUqVNpN0tKWWH/PxtjFllrkwqeW1HGbLsCa6216wCMMe8BQ4Eiw9Y1+9ZBjQTQU1IVypgxY1ixYgXp6emMHDlSQSsnraKEbSywOd/PW4BugSr8w8VbuHPqUq7tFU+bBtFUDQ+hakSI8+77c52q4YSkbIFX+kK7YXDeM+DxBqoJUsreeeed0m6ClHMVJWwL60YeNz5ijBkDjAFo3Lix34XHVo/EY+CtHzYUec4ZLevw1qgu0G0MzPkvZKTAsPHg9W/5NRGp2CpK2G4BGuX7OQ7YVvAka+14YDw4Y7b+Ft6taS1W/eMcDqdnczgj3ys9m5SMbG57dwnbD6Y7QwcDHobwaJj5CGSmwmUTITSypN9PRMq5ihK2C4DmxpgEYCtwOXBlICsI9XqoUSWMGlXCjvvss6Xb2LQv7Y8Dve+AiBj47M/wv4vhivcgIjqQzRGRcqZCTP2y1mYDtwDTgZXAVGvtr6XaqKRr4eLXYfN8mHg+pO4t1eaISOmqEGELYK39wlrbwlrbzFr7eGm3B4D2l8DwybD7N3jrHDh03MiGiFQSFSZsy6yWg2HE+3BoK7w5CPb+XtotqvR69ux5wnNGjx7NihXOzMF//vOfJ3191apVT61xZYw/a/r665577qFt27bcc889ASmvKMnJyX4tj1kcN/7+KsRDDacikA81jJm0kE370vjqjr5Fn7R1EUy+FKx1xnAbB2xmWrlSHh9qqFq1KocPH3b9moouOjqa3bt3H10/wi0TJkxg4cKFvPjii6dchr9/fyfzUIN6tsES2xmu/xoiq8OkC2DFJ6XdotL35f3w1rmBfX154ie183ot33//Pf379+eSSy6hVatWjBgx4uiiM/3792fhwoXcf//9HDlyhMTEREaMGHHM9YcPH2bAgAFH16/95JMT/50WtS5uXn0Ae/bsIT4+HoCcnBzuueceunTpQocOHXj11VcB2L59O3379iUxMZF27doxZ84ccnJyGDVqFO3ataN9+/Y888wzx9Wff/3d/N+lsPKAo2v65q33e8MNN9C2bVsGDhzIkSNHAFiwYAEdOnSgR48e3HPPPbRr1+64ei+44AJSU1Pp1q0bU6ZMKbIdxf2dLFiwgJ49e9KxY0e6du1KSkrKcfVkZmby8MMPM2XKFBITE5kyZUqR6wlPmDCBYcOGMXjwYJo3b8699957TFkPPfQQHTt2pHv37uzcubPYv1d/KGyDqVYzuH4m1O8AU0fCvJdPfI24asmSJTz77LOsWLGCdevW8cMPPxzz+ZNPPklkZCTJyclMnjz5mM8iIiL46KOPWLx4Md999x133XVXsSuE5V8X98MPP8Sf36zeeOMNYmJiWLBgAQsWLOC1115j/fr1vPPOOwwaNIjk5GSWLl1KYmIiycnJbN26leXLl7Ns2TKuvfZav/87FFZeQWvWrOHmm2/m119/pXr16kfX5L322mt55ZVXmDdvHl5v4Q/yTJs27eh/x+HDhxfblsL+TjIzMxk+fDjPPfccS5cuZebMmUeXvMwvLCyMRx99lOHDhx+tK2894SVLlvDoo4/y4IMPHj0/OTmZKVOmsGzZMqZMmcLmzc6zUampqXTv3p2lS5fSt29fXnvtNb//Wxalokz9KnU5uX4Ox1SpBSOnwYc3wPQH4MAmGPR45Xza7JwnS7sFdO3albi4OAASExPZsGEDvXv39utaay0PPvggs2fPxuPxsHXrVnbu3En9+vULPf9k1sXNM2PGDH755ZejvcCDBw+yZs0aunTpwnXXXUdWVhYXXnghiYmJNG3alHXr1nHrrbdy7rnnMnDgQL++B1BoeQUlJCQcPd65c2c2bNjAgQMHSElJOTqOfeWVV/LZZ5/5XW9hCvs7iYmJoUGDBnTp0gVwhiT8Vdx6wgMGDDi6qFCbNm3YuHEjjRo1Iiws7Oj+cJ07d+brr78u0XcC9WwDolWDaNbuPsy2A0f8uyA0Ei6dCN1vhvnjYOo1kJl24usk4PKPH+Zfg9YfkydPZvfu3SxatIjk5GTq1atX7NbWxfV6Q0JCyM3NBTimDGstL7zwAsnJySQnJ7N+/XoGDhxI3759mT17NrGxsVx99dVMmjSJGjVqsHTpUvr3789LL73E6NGji63HWktmZiZAoeUVVNh/q1O951NUO4qr51R3Zc5bT3j58uVHl8osri6A0NDQo/Wd7P8XRVHYBsClneMI8Rge+mgZmdm5/l3k8cLgf8LgJ+G3z31zcfe421A5JaGhoYXurnDw4EHq1q1LaGgo3333HRs3biy2nOLWxY2Pj2fRokUAx4xlDho0iHHjxh2tf/Xq1aSmprJx40bq1q3LDTfcwPXXX8/ixYvZs2cPubm5XHzxxTz22GMsXrz4uDbkr+eTTz45Wm5h5fmjRo0aVKtWjZ9++gmA9957z6/rimpHUVq1asW2bdtYsGABACkpKUUGYLVq1Y4Zzy0r6wkrbAOgUc0o/nZBW75btZub31lMVo6fgQvQ/UYY/jbsXA6vnQE7lrvXUDklY8aMoUOHDkdvkOUZMWIECxcuJCkpicmTJ9OqVatiyyluXdy7776bcePG0bNnz6ObTIIzBa1NmzZ06tSJdu3aMXbsWLKzs/n+++9JTEzk9NNP54MPPuD2229n69at9O/fn8TEREaNGsUTTzxxXBtuuOEGZs2aRdeuXZk/f/7R3SkKK89fb7zxBmPGjKFHjx5Ya/1a67eodhQlLCyMKVOmcOutt9KxY0fOPvvsIn+LOOOMM1ixYsXRG2RlZT1hTf0KoAk/rOdvn65g2OmxPHVpRzyek/i1Z+sieG8EpB9yFrBpfV5A21ZWlMepX4FUEdfFzftO4NxQ3L59O88991wptyo4NPWrlIzqlcBdZ7fgwyVbufXdJRzJPIl/RWM7ww3fQd1WMGUEzPqPMydXKpQxY8aQmJhIp06duPjii8t90AJ8/vnnx0wZ+8tf/lLaTSqT1LMNMGstr89Zzz+/XEnLetV48crTOa1uNf8LyEqHT2+DX6ZA24tg6MsQFhXwdpaWyt6zlcCZPn0699133zHHEhIS+Oijj4LWhsq4U0OZYYzhhr5NaVm/GndMSebc5+fywDmtuKZHvH/DCqERcNGrUK8tfP2I83jvFe9CTJz7jRcpRwYNGsSgQYNKuxl+0zCCS/q2qMNXt/ehR7Na/O3TFVz+2k+s35Pq38XGQK/b4cqpsH8DjD8DNs13tb3BVFl/m5KK5WT/P1bYuqhudARvjerCvy/pwMrthzj3+TlMXbjZ/7+kFgNh9EwIrwoTz4OFb5X7cdyIiAj27t2rwJVyzVrL3r17iYiI8PsajdkGyfaDR/jzlGR+WrePM1vV5cEhrfwfy03bBx+Mht+/gcQRcO5/y+3uD1lZWWzZsqXYyf8i5UFERARxcXGEhh679VVRY7YK2yDKybW89cN6npu5hrSsHK7s2pi7BragetTxuz8cJzcHZv3LedVrD8MnQc2m7jdaRE6KwraA0gjbPHsPZ/DcN2uYPH8TNaJC+et5bbigY0P/HkdcPcNZV8FaGPYqtDzH/QaLiN80z7YMqVU1nEeHtmPaLb2IrRHF7e8lc/n4n0je7McizS0GwthZUDMe3r0cvnnU6fWKSJmmsC1FbRvG8OGNPXnswnas3XWYC1/6gZsmL2Ld7hMsWlwjHq6bAadf7Wyb/r9hWldBpIzTMEIZcTgjm9dmr+P1OetIz87l0s5x3DagOQ2rn+BG2OJJ8PndUKW2s5JYoy7BabCIFEpjtgWUtbDNs+dwBi99t5bJP20CA1d3b8JN/ZtRq2oxW4lsW+Is03hoG5z5V+h5G3j0S4tIaVDYFlBWwzbPlv1pPP/NGt5ftIXIUC/X905gdN+mREeEFn7BkQMw7VZYOQ2aDXCeQqtaJ7iNFhGFbUFlPWzzrN11mGe+Xs3ny7YTExnKjf2bMbJHPJFhhezsYC0sfAO+etDZ62zYa9C0X/AbLVKJKWwLKC9hm2f51oM8NWMV36/aTZ1q4Yzt25QruzUmKqyQ5S12LIP/uxb2roW+90C/+8CrZTBEgkFhW0B5C9s8Czbs4+kZq5m3bi81q4Rxfe8Eru7R5PjhhYzD8MU9sPQdaNwTLn4dYmJLp9EilYjCtoDyGrZ5Fm3cx4vfruW7VbupFhHCyB7xXNc7gZpVCjyNtvQ9+OxOCAmHC8dBy8Gl02CRSkJhW0B5D9s8y7ce5KXv1vLl8h1EhnoZ0a0xY/o2pW50vgUy9qxxhhV2LoOuY+DsR8vt2goiZZ3CtoCKErZ51uxM4eXvf+eT5K2EeD1clhTH2L7NaFTTt/B4VjrM/Juzm2+dVs7NswYdSrXNIhWRwraAiha2eTbuTeWVWb/z/qItWAsXnh7Ljf2b0ayOs0cUa2fCxzfBkf3OnNwet2hOrkgAKWwLqKhhm2fbgSOMn72Od3/eRGZOLme3rsfYfk3p3KQmpO51tt757TNI6AcXvQLRDUu7ySIVgsK2gIoetnn2HM5g4o8bmDRvIwePZNG5SQ3G9G3K2a3q4kl+G766H7xhcP5z0PbC0m6uSLmnsC2gsoRtnrTMbKYu2Mzrc9ezZf8Rmtapwg19mjKs8RHCp411HvlNvArOeRLCT2KDShE5hsK2gMoWtnmyc3L5YvkOxs/+neVbD1G7ajjXdovlupwpRM5/Dqo3dh71bdy9tJsqUi4pbAuorGGbx1rLvN/38ursdcxavZuIUA93tdzLqJ1PEpqyBXrdBv0fdHb7FRG/KWwLqOxhm9/qnSm8PmcdHy/ZRlhuKi/X/pC+KZ9j67TGDHsVGnQs7SaKlBsK2wIUtsfblZLOpB838r/5G0lM/5mnI96guj2E7Xcf3j53an0FET8obAtQ2BYtLTObDxZtYeqcXxidMo6h3h/ZFd2WqOGvUzW2TWk3T6RM0x5k4reosBCu7hHPx3efT8Tlb/FMzIOEHtxIyPi+TH/9YTbvPcG2PSJyHPVsxS8rVq/CfnIbbVN/Yl5uGz6L/wsXntmDpCY1/NsVWKSS0DBCAQrbU2AtB358k8hv/0JOTg7/zLqS5fWHcW2fZgxp34BQr35RElHYFqCwLYEDm8n55Ba8679nsbc9t6VdT050Y67pEc+VXRsTE1XE1j0ilYDCtgCFbQlZC4snYac/RE5uDpOqXMdjO7sTGRbKpZ3juLZXAvG1q5R2K0WCTmFbgMI2QA5shk9vh9+/IbVBD56tcjsTVlqycy1nt67H6D5N6RKvcV2pPBS2BShsA8haWPI2TH8IcrNJ6f0XXj1yJpN/3sz+tCzax8Ywuk+CxnWlUlDYFqCwdcHBLU4vd+1MaNKb9CHP8cGGUN6Yu551u1OpHx3BqF7xXNFF47pScSlsC1DYusRaSJ4MXz0Audlw5l/J7TKG79fu5fU56/nx971EhXk1risVlsK2AIWtyw5uhc/+DGumQ1xXGPoi1GnJim2HeGPueqYt3apxXamQFLYFKGyDwFpY9n/w5b2QmQr97oVed4A3lF2H0pk0z1mH4UBaFh3iYhjdpynntKuvcV0p1xS2BShsg+jwLvjiHljxMdRvD0NfOrqS2JHMHD5YvIU3565n3Z5UGsY447rDuzQmJlLjulL+KGwLUNiWgpWfwud3Qeoe6HU79Lvv6Hq5ubmW71bt4vU565m3bi9Vwrxc1qUR1/VK+GOHYJFyQGFbgMK2lBzZ70wRS54MtVvABS9C427HnLJ860HenLueaUu3kWstg9rWZ3SfBGezSpEyTmFbgMK2lK2dCZ/e4UwX6zbW2VY9vOoxp+w4mM7EeRuY/NNGDqVnk9ioOqP7JDC4bX1CNK4rZZTCtgCFbRmQkQIz/w4LXoOYxnD+M3DaWcedlpqRzQeLt/DG3PVs3JtGXI1IbujTlMuSGhEZ5i2FhosUTWFbgMK2DNn0E0y7Ffashg6Xw+AnIOr4IYOcXMvMlTt5ddbvLN50gBpRoYzsGc81PeKpWSWsFBoucjyFbQEK2zImKx1m/wd+eBYiqsOQf0PbYVDE3NsFG/bx6qzfmblyFxGhHoYnNWJ0n6a6mSalTmFbgMK2jNqxHKbdAtuWQItz4Nz/Qkxskaev2ZnC+Nnr+Dh5Kzm5lnM7NGRs36a0i40JYqNF/qCwLUBhW4blZMP8cfDt4+AJgbP/Dp2vBU/RN8V2HEznrR/WM3n+Jg5nZNOneW3G9m1Gr9Nq6ck0CaoytQeZMeZSY8yvxphcY0xSgc8eMMasNcasMsYMynd8sO/YWmPM/fmOJxhj5htj1hhjphhjNHhX3nlDoOetcNOPEHs6fH4nTDwP9qwt8pL6MRE8MKQ1P9x/JvcNbsVvO1K46o35nPfCXKYt3UZ2Tm4Qv4DI8UqlZ2uMaQ3kAq8Cd1trF/qOtwHeBboCDYGZQAvfZauBs4EtwALgCmvtCmPMVOBDa+17xphXgKXW2nEnaoN6tuWEtbDkfzDjIWdct/990PM28Bb/dFlGdg4fL9nKq7PXsW53qmYwSNCUqZ6ttXaltXZVIR8NBd6z1mZYa9cDa3GCtyuw1lq7zlqbCbwHDDXO74dnAu/7rp8IXOj+N5CgMQY6XQ03/wwtBsE3j8L4M2Dr4mIvCw/xMrxLY2b+uR+vXt2ZOtXCeWTar/T617e8+O0aDh7JCtIXEHGUtZnhscDmfD9v8R0r6ngt4IC1NrvAcaloqtWH4W/D8P9B6m54fQDM+AtkphV7mcdjGNS2Ph/e2JOpY3vQMS6Gp2aspteT3/LElyvZlZIepC8glV2IWwUbY2YC9Qv56CFr7SdFXVbIMUvh/yjYYs4vqk1jgDEAjRs3Luo0Kctanw/xfeDrh+HHF5z1Fs5/Hpr2K/YyYwxdE2rSNaErK7YdYtys33lt9jre+mEDl3aOY2zfZjSupWlj4h7XerbW2rOste0KeRUVtOD0TBvl+zkO2FbM8T1AdWNMSIHjRbVpvLU2yVqbVKdOnVP5WlIWRFaHC56HkZ+C8cCkC+CTW5x1F/zQpmE0L1xxOt/e1Z+LO8Xyfwu3cMZ/v+f295bw245DLjdeKquyNowwDbjcGBNujEkAmgM/49wQa+6beRAGXA5Ms87dve+AS3zXjwSKC3OpSBL6wo0/OmvkJr8DL3WDFf7/9cfXrsITwzow574zuK5XPF+v2MngZ+cweuICFm30L7hF/FVasxEuAl4A6gAHgGRr7SDfZw8B1wHZwB3W2i99x4cAzwJe4E1r7eO+401xbpjVBJYAV1lrM07UBs1GqGC2JTuP/O74BVqd5zwMUa2wUayi7U/NZOK8DUz4cQMH0rLollCTm884jT7Na2uurvhNDzUUoLCtgHKyYd6L8P0T4A2HQf+A068u8pHfoqRmZPPuz5t4fc56dhxKp11sNDf1P41Bbevj9Sh0pXgK2wIUthXY3t9h2m2wca4z1HD+c1Cz6UkXkzdX95VZ61i/J5Wmtavwp37NuPD0WMJCytoInJQVCtsCFLYVXG4uLJ4AMx727fL7EHS/CTwn/0BDTq7ly+Xbefm731mx/RBxNSK5fUBzLjo9VuvqynEUtgUobCuJg1udx31XfwUNOzm7/NZre0pFWWv5ftVunv56Ncu2HqRp7SrccXYLzmvfAI+GF8RHYVuAwrYSsRaWfwBf3gfpB6D3ndD3bggJP8XiLDNW7OTpGatZtTOFVvWr8eezWzCwTT3dSBOFbUEK20oodS9MfwB+mQK1Wzq93EZdT7m43FzLZ8u28+zXq1m3J5UOcTHceXYL+rWoo9CtxBS2BShsK7E1Xzv7nx3a6ozjnvkXCDv1p8eyc3L5aMlWnvtmDVv2HyGpSQ3uGtiSHs1qBbDRUl4obAtQ2FZy6Ydg5t9g4RtQIwEueAES+pSoyMzsXKYu3MwL365h56EMep9WmzsHtqBT4xqBabOUCwrbAhS2AsD6Oc7DEPvXQ9L1zkLl4dVKVGR6Vg6T52/i5e/Wsjc1kwGt6vLns1to94hKQmFbgMJWjspMg+8eh3kvQUwcnP9sobv8nqzUjGwmztvAq7PWcfBIFkPa1+fPZ7Wgeb2ShbmUbQrbAhS2cpzNPzsL2uxZBYkjYNDjEFnyIYBD6Vm8MWc9b8xdT2pmNhcmxnLHWc1pUqtKABotZY3CtgCFrRQqKx1m/xvmPgtVasN5z0CrcwNS9P7UTF6dvY4JP64nK8dyVbfG3HFWC2poG/YKRWFbgMJWirUt2enl7lwG7S+Fc/4NUTUDUvSulHSe/2YN78zfRNXwEG4b0JxresTrEeAKQmFbgMJWTignC+b8F2b/B6JqwXnPQqshASt+9c4UHv98JbNW7ya+VhQPDGmtByMqgDK1B5lIueANhf73ww3fQZU68N4V8OFYvxcpP5EW9aox8bquTLi2C6FeD2PfXsQVr/3E8q0HA1K+lC3q2Yr4IzsT5jwFs5+CqnWdlcRaDApc8Tm5vLdgM09/vZr9aZlc3CmOewa1pF50RMDqkODQMEIBCls5JduS4eObYNevvhkL/3S26QmQQ+lZvPTtWt76YQNej+FP/Zoxpm9Tbb9ejihsC1DYyinLzoBZ/4a5z0DVes7TZ81LPi83v01703jyq5V8sWwH9aMjuHdwSy5MjNXqYuWAxmxFAiUkHAb8FUbPhIhomHyx8xRaRkrAqmhcK4qXR3Tm//7Ug7rR4dw5dSkXvvwDP6/fF7A6JLjUsxUpiewMZxueH56DmEZw0SvQpGdAq8jNtXycvJV/f7WKHYfSOaddfR44p7W2Xi+jNIxQgMJWAmrTT/DRWNi/EXrdBmc8dMrr5RblSGYO42ev45VZv5OTa7m2Vzw3n3ka0RGhAa1HSkZhW4DCVgIu4zDMeAgWTYC6bWHYeKjfLuDV7DyUzn+mr+KDxVuoERXGA+e04pLOcZqfW0ZozFbEbeFVnSlhV06F1N0wvr9zEy03J6DV1IuO4KlLO/LpLb1pWrsK97z/C9e8+TOb96UFtB4JLPVsRdyQuhc+uwNWToNG3Z2x3JoJAa8mN9cyef5GnvzyN3It3DOoJSN7xmvL9VKknq1IMFWpBZdNgovGw66VMK4XLJro7IcWQB6P4eoe8cy4sx/dm9bk0c9WcPG4H1m9M3AzIyQwFLYibjEGOg6Hm36EuM7w6W0w5SpIC/z0rdjqkbw5qgvPXZ7Ixr2pnPv8HJ6duZrM7NyA1yWnRmEr4raYOLj6Exj4D1g93enlrp8d8GqMMQxNjGXmnf0Y0r4Bz85cw3kvzGHJpsCs5SAlo7AVCQaPB3re6jwIERYFEy9w9kDLyQp4VbWqhvPc5afz5qgkUtKzGTbuRx79dAVpmdkBr0v8p7AVCaaGiTB2NnS62pmp8MZA2Pu7K1Wd2aoeM/7clxHdGvPmD+sZ+Mxs5q7Z40pdcmIKW5FgC6virKdw2STYtw5e6QNLJgf85hlAtYhQ/nFhe6aO7UGY18NVb8znnv9bysG0wPeopXgKW5HS0mYo3PgDxHaCT26C96+FIwdcqaprQk2+uL0PN/VvxodLtjLg6Vl8uWy7K3VJ4RS2IqUpJg6u+QQGPAwrP3V6uVvcmf8dEerl3sGtmHZLL+rHhHPj5MWMfXshuw6lu1KfHEthK1LaPF7ocxdcNx0M8OYgZ1t1lx44atswho9v6sX957Ti+1W7GfD0LKYs2ERlfcApWBS2It2x2I0AACAASURBVGVFXJJz86zFYJj+ILx3ZcC24CkoxOvhT/2a8dUdfWnTIJr7PljGre8u4XCGZiy4RWErUpZE1oDh/4PBT8Kar+GVvq4NKwAk1K7Cuzd0597BLfli2XYueHEuq3bo6TM3KGxFyhpjoPuNzrACwJuDYd7Lrg0reDyGm/qfxuTR3Tl0JJuhL83lw8VbXKmrMlPYipRVcZ3hT7Oh+UCY/gC8N8K1YQWAHs1q8cVtvekYV507py7lgQ9/IT0rsCuWVWYKW5GyLLIGXD4ZBj0Ba6Y7wwpbF7tWXd3oCCaP7sZN/Zvx7s+buXjcj2zcm+pafZWJwlakrDMGetzkG1awzrDC4rddqy7E6+Hewa14c1QSW/Yf4bwX5jL91x2u1VdZKGxFyou4JBgzCxp3h2m3wGd/huxM16o7s1U9PrvVWaB87NuLePzzFWTlaBWxU6WwFSlPqtSCqz6EXnfAwjdhwhA4tM216hrVjGLqn3pwTY8mvDZnPVeM/4kdB/UQxKlQ2IqUN94QOPvvcOlE2LkCXu0HG390rbrwEC+PDm3H81eczorthzj3+Tla0OYUKGxFyqu2F8IN30JENEw8H356xbXpYQAXdGzItFt6U6tqGFe/OZ/nv1lDbq6eOvOXwlakPKvbygnc5gPhq/uc7dQz3dv48bS6Vfn45l5clBjL01+vZtSEBexLdW/cuCJR2IqUdxExMHwynPEX+GWqbxzXvRW9osJC+O9lHXliWHt+WreXc5+fw6KN2g3iRBS2IhWBxwP97oEr3oXdq+G1M2BbsmvVGWO4omtjPryxJ6FeD8NfncfUhZtdq68iUNiKVCQtz4HrZ4AnxJmPu2Kaq9W1i43h01t70/O02tz7/i9M+GG9q/WVZwpbkYqmfjtnHLd+O5h6Ncx+ytUbZzGRobx2TWcGta3H3z5dwbjv3dnmp7xT2IpURFXrwsjPoP2l8O1jzo2zLPfmx4aHeHnxyk4MTWzIv776jadnrNL6uAWElHYDRMQloREw7DWo3RK++wfs3+DcSKtax53qvB6eviyRyFAvz3+7lrTMHB46tzXGGFfqK2/UsxWpyIxxbpxdOgG2/wKvnencQHOJ12P450XtGdUzntfnrucvHy/XXFwfha1IZdD2Irj2C8g+Am8OhM0/u1aVx2N45Pw23Ni/GZPnb+Lu95eSrTUVFLYilUZsJ2emQmQNmHgB/PaFa1UZY7h3UEvuOrsFHy7eyu3vJZOZXbkDV2ErUpnUbArXzYC6rWHKCFg0wbWqjDHcOqA5fzm3NZ8v286N/1tUqRcjV9iKVDZV68DIT6HZAPj0dvjuCVenho3u05THLmzHN7/tYvTEhaRlVs5NJRW2IpVReFXnabPEETDrSSd0c9wLwau7N+GpSzvy4+97GPnmz6SkZ7lWV1mlsBWprLyhMPQl6HM3LJ4IU66CrCOuVXdJ5zheuKITSzYdYMTr8zmQVrkWsFHYilRmxsCAv8KQp2D1VzD5Usg47Fp153ZowCtXdea37SlcPv4ndqdkuFZXWVMqYWuM+Y8x5jdjzC/GmI+MMdXzffaAMWatMWaVMWZQvuODfcfWGmPuz3c8wRgz3xizxhgzxRgTFuzvI1Ludb0BLnoVNv4A/xsG6Qddq+qsNvV4Y1QSG/amMnz8vEqz80Np9Wy/BtpZazsAq4EHAIwxbYDLgbbAYOBlY4zXGOMFXgLOAdoAV/jOBfgX8Iy1tjmwH7g+qN9EpKLoOBwueQu2LnKmhqXtc62qPs3rMOm6buw6lMGot34mNaPi3zQrlbC11s6w1ub91/0JiPP9eSjwnrU2w1q7HlgLdPW91lpr11lrM4H3gKHGeQ7wTOB93/UTgQuD9T1EKpy2FzqP9O5aCRPOg8O7Xauqa0JNXhrRidU7U7hr6tIK/6RZWRizvQ740vfnWCD/ophbfMeKOl4LOJAvuPOOF8oYM8YYs9AYs3D3bvf+JxIp11oOhiunwL518NY5rm4o2a9FHR4c0pqvft3Bc9+sca2essC1sDXGzDTGLC/kNTTfOQ8B2cDkvEOFFGVP4XihrLXjrbVJ1tqkOnXcWYxDpEJodgZc/SGkbHcCd/9G16q6vncCl3SO47lv1vD5L+7tMFHaXFv1y1p7VnGfG2NGAucBA+wfa7FtARrlOy0OyPtntbDje4DqxpgQX+82//kiUhJNesI1nzg3zCacC9d+CdUbnfi6k2SM4fGL2rFu92Hu+r9kmtSKol1sTMDrKW2lNRthMHAfcIG1Nv/udNOAy40x4caYBKA58DOwAGjum3kQhnMTbZovpL8DLvFdPxL4JFjfQ6TCi0tyAjf9kLODr0t7m4WHeHnl6s7UiApjzKSFFXJKWGmN2b4IVAO+NsYkG2NeAbDW/gpMBVYAXwE3W2tzfL3WW4DpwEpgqu9ccEL7TmPMWpwx3DeC+1VEKriGp8NV78PhXTBpKKTucaWautUieO2aJPalZfKn/y0iI7tiraNgKutq6klJSXbhwoWl3QyR8mPDXPjfxVCrOYycBlE1Xanms1+2ccs7S7gsKY5/Xdyh3C0+boxZZK1NKni8LMxGEJHyIL43XP4O7FnlhG76IVeqOa9DQ2478zSmLtzCWz9scKWO0qCwFRH/nTYALpsEO36Bdy6DzFRXqrnjrBYMaluPf3y+gtmrK8Y0TYWtiJyclufAxa/D5vnw7hWQHfibWR6P4enLEmlRrxq3vLOYdbvdW68hWBS2InLy2l7krBi2fhZ8fBPkBn4XhirhIbx2TRIhXg+jJy3k4JHyvSyjwlZETk3ilTDgEVj+Psx8xJUqGtWMYtyITmzam8Zt7y4hpxw/0quwFZFT1/vP0GU0/Pg8zH/VlSq6Na3FYxe2Y9bq3Tz55UpX6ggG154gE5FKwBg459+QsgO+vA+q1Yc2Q0983Um6omtjftt+iNfmrKdzk5oMblc/4HW4TT1bESkZj9e5YRbXBT64ATbOc6Wav57Xhlb1q/HYZyvK5caRClsRKbnQSGelsOqN4N3LYfeqgFcR4vXwyPlt2XrgCK/OWhfw8t2msBWRwIiqCVd9AN4weGe4K4uP92hWiyHt6zNu1lq2HXBvvzQ3KGxFJHBqxMPlk+HQVnj/Old27H1wSGushSe+/C3gZbtJYSsigdWoK5z3DKz7zpUpYXE1ohjbrxmfLt3Gz+vd27on0BS2IhJ4p18F3f4E816E5HcDXvyf+jWlQUwEf5v2a7mZe6uwFRF3DPwHJPSFT2+HLYsCWnRUWAgPDGnNiu2HmLpw84kvKAMUtiLiDm8oXDrRmXs7ZYQzFzeAzu/QgK7xNfnP9FXl4lFeha2IuCeqJlzxrrMc49RrICdwoWiM4eHz27A/LZPnZpb9zSIVtiLirnptYegLziph3z0e0KLbxcZweZdGTJq3gbW7UgJadqApbEXEfe0uhs7XwtxnYO3MgBZ998CWRIZ5efSzlZTlnWcUtiISHIOfgLpt4cOxAd04slbVcO44qwWzV+/mm5W7AlZuoClsRSQ4QiPh0gmQlQYf3gC5gVvf4JoeTWhWpwqPfb6izG4UqbAVkeCp0wLOfRo2zIFZ/w5YsaFeDw+f35aNe9PK7L5lClsRCa7EK6DjlTDrX7B+TsCK7deiDme1rssL36zhQFpmwMoNFIWtiATfkP9ArWbw8Y0B3aX3roEtSc3MYcqCsvegg8JWRIIvvCpc+IqzYM2MhwJWbOsG0XRvWpNJ8zaSnRP4fdFKQmErIqWjURfoeRssngRrvg5YsaN6JrD1wBFmlrGZCQpbESk9/R+AOq1g2q1wZH9AijyrdV1iq0cy4cf1ASkvUBS2IlJ6QiPgolfg8C746oGAFBni9XBNjyb8tG4fK7cHbjy4pBS2IlK6Gp4Ofe6Cpe/Cb18EpMjhXRoREephQhmaBqawFZHS1/ceqNcePrsD0g+WuLjqUWEM6xTHx8lb2ZdaNqaBKWxFpPSFhMEFzznDCd89EZAiR/WMJyM7l/cWbApIeSWlsBWRsiG2MyRdBz+/CtuXlri4FvWq0eu0WrxdRqaBKWxFpOwY8FeIqgWf3Qm5JQ/IUT0T2H4wnRkrdgagcSWjsBWRsiOyhrOdztaFsHhiiYs7s1VdGtWMLBM3yhS2IlK2dBgOTXrDzL9B6p4SFeX1GEb2iOfnDftYvrXkN95KQmErImWLMXDufyHzcEC2Qr80qRGRoV4m/rih5G0rAYWtiJQ9dVs5W6EvmQw7lpWoqJjIUC7uHMsnS7ex93BGgBp48hS2IlI29b0bImLg64dLXNTIHvFkZufyzvzSmwZWZNgaY37x4/VNMBsrIpVIZA3ody/8/m2J9y1rXq8a3ZvWZNrSbQFq3MkLKeYzLzCkmM8NMC2wzRERyafLaPh5PMx4GJqeAR7vKRd1Zqu6/POL39h+8AgNYiID2Ej/FDeMMNZau7GY1wbgpiC1U0Qqo5BwGPAI7PrVWTuhBPo0rwPAnDUlm+FwqooMW2vt3BNd7M85IiIl0vYiiE2Cb/8BmWmnXEyr+tWoUy287IVtHmPMecaYJcaYfcaYQ8aYFGNM2Vm3TEQqNmOcBx1StsPCN0pQjKFP89rMXbObnFwbwAb6x5/ZCM8CI4Fa1tpoa201a220y+0SEflDkx7QtD/88FyJerd9m9dhf1oWv24L/gMO/oTtZmC5tTb4/xSIiOTpdz+k7oaFb55yEb2b1wZg9urdgWqV3/wJ23uBL4wxDxhj7sx7ud0wEZFjNOkBCf1K1LutXTWctg2jmV0K47b+hO3jQBoQAVTL9xIRCa7+90PqLlj01ikX0bdFHRZv3M/hjOwANuzEiptnm6emtXag6y0RETmRJj0hvo/Tu026DkJPfr5sn+a1Gff978z7fS9nt6nnQiML50/PdqYxRmErImVD//vh8E5Y/PYpXd65SQ0iQ73MWRPccVt/wvZm4CtjzBFN/RKRUhffG+K6wE8vQ27OSV8eHuKlR7NaQb9JdsKw9U318lhrIzX1S0TKhB43w/71sOrUduPt07w2G/amsWnvqU8jO1nFLURT/0QX+3OOiEjAtTofqjeGeS+d0uV9W/ge3V0bvN5tcT1bf/7JCMwm7yIiJ8MbAt1uhE3zYMuik768ae0qxFaPDOpQQnFh29E3RlvUKwUI3q08EZH8Tr8KwqPhp5Pv3Rpj6NuiNj+u3Ru0nXeLW4jG6xujLepVzVobG5RWiogUFBENna6BXz+GA5tP+vLuTWuRkpHNml2HXWjc8bRTg4iUX93Ggs2FxZNO+tK2DZ37/L9uC87kKoWtiJRf1RtD87NhyduQc3JPhCXUrkpkqDdoi9IobEWkfOs8yll+cc30k7rM6zG0alCtbPVsjTFeY0xDY0zjvFdJKjXGPObbwyzZGDPDGNPQd9wYY543xqz1fd4p3zUjjTFrfK+R+Y53NsYs813zvDHGlKRtIlLONB8E1RrAogknfWnbhtGs3HaI3CCsb+vP4uG3AjuBr4HPfa/PSljvf6y1Hay1ib6y8rbPPAdo7nuNAcb52lATeAToBnQFHjHG1PBdM853bt51g0vYNhEpT7whcPrVsOZrOHByu+e2bRhDSkY2m/e7/3CDPz3b24GW1tq21tr2vleHklRqrc3fb68C5P2zMhSYZB0/AdWNMQ2AQcDX1tp91tr9OME/2PdZtLV2nm+93UnAhSVpm4iUQ52udt5Pcr2EYN4k83fx8ICPIBtjHjfGbAZG8EfPNtZXX54tvmPFHd9SyPGi6hxjjFlojFm4e3fwFw8WEZdUbwynnQXJkyHX/3mzLepVw+sxQblJVtzjunmLhK8Dvj/ZxcONMTONMcsLeQ0FsNY+ZK1tBEwGbsm7rJCi7CkcL5S1dry1Nslam1SnTp0TfQURKU86Xg6HtsKmH/2+JCLUS/O6VVkRhJ5tcevZ5i0Qvsn3CvO9oJhAy2OtPcvPNryDMw78CE7PtFG+z+KAbb7j/Qsc/953PK6Q80Wksmk5BMKqwi9TnJXB/NSmYTRzg7BzQ3FPkP3dWvt3YEXen/MdW1mSSo0xzfP9eAHwm+/P04BrfLMSugMHrbXbgenAQGNMDd+NsYHAdN9nKcaY7r5ZCNcAn5SkbSJSToVFQavz4NdPICvd78vaNoxhV0oGu1MyXGycf2O2D/h57GQ86RtS+AUnOG/3Hf8CZ9hiLfAacBOAtXYf8BiwwPd61HcM4Ebgdd81vwNflrBtIlJedbgMMg7Cmhl+X/LHTTJ3x22LHEYwxpwDDAFijTHP5/soGijR5j3W2ouLOG5xFisv7LM3geO21bTWLgTalaQ9IlJBJPSDKnWdoYQ2F/h1SesGf8xI6N+yrmtNK65nuw1YCKQDi/K9puFMxRIRKVu8IdD+Eqdne+SAX5fERIbSqGak6zfJihuzXWqtnQicBrwLLAEWA5/55rqKiJQ9bS+CnMyTG0poEOP6MII/Y7Zn44yFPg+8CKz1DTGIiJQ9sUlQtT6s/NTvS9o2jGbD3jRS0rNca5Y/Yfs0cIa1tr+1th9wBvCMay0SESkJjwdanQtrZ0LWEb8uaVHfmem6fk+qe83y45xd1tq1+X5eB+xyqT0iIiXX+nzISoPfv/Xr9IYxkQDsOOj/lLGT5U/Y/mqM+cIYM8q32tanwAJjzDBjzDDXWiYicqrie0NEdb+HEurHRACw3cWwLe4JsjwROKt+9fP9vBuoCZyP8yTZh+40TUTkFHlDoeU5sOpLyMlyfi5GrSphhHpN6YattfZa12oXEXFLyyGw9F3YPP+Ej+96PIb6MRHsOOjfGO+p8Gc92xbGmG+MMct9P3cwxvzFtRaJiARC0/7gCXFulPmhQXSkqz1bf8ZsX8N5PDcLwFr7C3C5ay0SEQmEiGho1B3W+Be29WMiSj1so6y1Pxc4VqLHdUVEgqL5WbBzGRzafsJTG8REsONgOs6qAYHnT9juMcY0w7esojHmEuDELRcRKW2n+VZ6/f2bE57aICaCzJxc9qVmutIUf8L2ZuBVoJUxZitwB85KWyIiZVu9ds5mkGu+PuGp9X1zbd0aSvBnNsI64CxjTBXAY61NcaUlIiKBZgycNsCZb5ubAx5vkac2yDfXtl1sTMCbUtwSi4VufZO3U7i19umAt0ZEJNAS+sGS/8GOZdAwscjTGlR3wtat6V/FDSNU872ScIYN8jZY/BPQxpXWiIgEWt4c2w1zij2tdpVwQjzuPdjgz7Y4tYFO1tq7rLV3AZ05dt8vEZGyK7oh1GwGG+YWe5rHY6gXHeHa+gj+3CBrDOS/PZcJxLvSGhERNyT0gY0/Qk7xs1YbxESwrRSGEfK8DfxsjPmbMeYRYD4w0ZXWiIi4Ib4PZByCHUuLPa1B9cjS69laax8HrgX2AweAa621T7jSGhERN+SN264vfty2fnQ4Ow6V0tQvAGvtYpwtcUREyp9q9Z1x280FH4Y9Vs0q4aRn5XIkM4fIsKKniZ0Kf4YRRETKv0ZdnRXAinkct2YVZynGfWmBf4pMYSsilUNcF0jbA/vXF3lKjagwAPa78MiuwlZEKodG3Zz3zQuKPKVmFSds3VgfQWErIpVD3dYQVs0ZSihCDV/Y7tcwgojIKfJ4Ia4zbCn6JlnNKPVsRURKLq4L7PwVMtMK/Tg6MhSP0ZitiEjJNEgEmws7lxf6sddjqB4VptkIIiIlkrfq17bkIk+pERXK/tSsgFetsBWRyiM6FqJqw/aiH9utWSWMvakZAa9aYSsilYcxTu92e3E92zD1bEVESqxBIuxaCVmFr+5Vs4rGbEVESq5hItgc2Lmi0I9rVAljf2pmwHfZVdiKSOVSr53zXsSMhJpRYWTnWlIyil/79mQpbEWkcqneBEKrwK6ie7YQ+Lm2ClsRqVw8HufR3Z2/Fvrx0ZW/FLYiIiVUr40TtoWMy0aFOct8H8nKCWiVClsRqXzqtoUj++DwzuM+Cg9xYjEjOzegVSpsRaTyqdfGeS9kKCE8xNmhISNLYSsiUjJ1fWG7e9VxH4WH5vVsNYwgIlIyUbUgsgbsWX3cRxpGEBEJFGOgdgvYu/a4j44OIyhsRUQCoFbzwnu2ecMImo0gIhIAtZs7sxGOHDjmsIYRREQCqXYL573AUEKYV2ErIhI4tZs773vWHHPYGEN4iEezEUREAqJ6EzAe2L/+uI/CQzyaZysiEhAhYRATB/vWHfdReKhXwwgiIgFTIwH2FdGz1TCCiEiA1EwoehhBPVsRkQCp2RTS9kL6wWMOh4d4NWYrIhIwNRKc9wJDCeGhGkYQEQmc6o2d94Objzms2QgiIoGUF7YHNh1zODzES0aOwlZEJDAia0BYVThQWM9WwwgiIoFhDMQ0On4YIdRLpmYjiIgEUPVGcGDjMYcq3NQvY8zdxhhrjKnt+9kYY543xqw1xvxijOmU79yRxpg1vtfIfMc7G2OW+a553hhjSuO7iEg5FdOo8GGEijIbwRjTCDgbyD8yfQ7Q3PcaA4zznVsTeAToBnQFHjHG1PBdM853bt51g4PRfhGpIGLiIP0AZBw+eqiizbN9BrgXyL+X8FBgknX8BFQ3xjQABgFfW2v3WWv3A18Dg32fRVtr51lrLTAJuDC4X0NEyrXoWOc9ZfvRQ8482woQtsaYC4Ct1tqlBT6KBfL357f4jhV3fEshx4uqd4wxZqExZuHu3btL8A1EpMKIbui8H9p69FB4iIfMnFxyc20RF528kICVVIAxZiZQv5CPHgIeBAYWdlkhx+wpHC+UtXY8MB4gKSkpcP8VRaT8Ohq2+Xq2vn3IMnNyifB4A1KNa2FrrT2rsOPGmPZAArDUdy8rDlhsjOmK0zNtlO/0OGCb73j/Ase/9x2PK+R8ERH/VGvgvBfo2QJkZOUSERqYsA36MIK1dpm1tq61Nt5aG48TmJ2stTuAacA1vlkJ3YGD1trtwHRgoDGmhu/G2EBguu+zFGNMd98shGuAT4L9nUSkHAuLch5uKDBmCwR0RoJrPdtT9AUwBFgLpAHXAlhr9xljHgMW+M571Fq7z/fnG4EJQCTwpe8lIuK/ag3h0B+/FLuxnXmph62vd5v3ZwvcXMR5bwJvFnJ8IdDOrfaJSCVQrZ6z067PHzvsBq5nqyfIRESq1oPDu47+mBe26QGca6uwFRGpWtfp2VpnklJ4aOCHERS2IiJV60NOJhzZD0CY14nGQC5Go7AVEala13n3DSWEep0p/Nm5ClsRkcCpWs95P7wDAK8nL2wD9+yTwlZEJK9nm7oHgBCPE405OQpbEZHAiartvPvCVj1bERE3RNYA44E0X8/WN2abo7AVEQkgjweiakGqsxrgHz1b3SATEQmsqNr5xmzVsxURcUeV2hqzFRFxXZXaf4zZ5s1GUNiKiARYZE1IcxYTVM9WRMQtUTWdjR9zc/8Ys83RDTIRkcCKrAk2F9IP4PWqZysi4o6oms77kf2ajSAi4ppIX9im7dOYrYiIa472bPdpNoKIiGsiqjvv6QfxdWzJ1g0yEZEAi4hx3o8cwBhDiMeQY9WzFREJrMi8nu0BwJlrqzFbEZFA84ZCaBU44oRtiMdoPVsREVdExED6QUA9WxER90RWPzqMEOL1aDaCiIgr1LMVEQmC8GpHwzbEY8jR4uEiIi4Ij4aMFEA9WxER90REQ8YhIK9nq7AVEQm88GhIPwTWqmcrIuKa8GqQmwXZGYR4PJpnKyLiirxHdjMOqWcrIuKasKrOe0YKIV7NRhARcUe4L2wzD6tnKyLimryebWaqZiOIiLgmvJrznqGerYiIe472bFPwqmcrIuKSvDHbjMMYTECLVtiKiOQJq+K8Z6YGvGiFrYhInlBf2GYpbEVE3BMSBp5Q9WxFRFwXFgWZaQEvVmErIpJfaBUNI4iIuC6sioYRRERcp2EEEZEgCI2C7CMBL1ZhKyKSX2ikerYiIq4LjYIs9WxFRNwVGglZ6tmKiLgrNFI9WxER12kYQUQkCEIiNBtBRMR1oZGQk4mxOQEtVmErIpJfSDgAYTYzoMUqbEVE8guJBCBUYSsi4qLQCADCbEZAi1XYiojk5+vZVohhBGPM34wxW40xyb7XkHyfPWCMWWuMWWWMGZTv+GDfsbXGmPvzHU8wxsw3xqwxxkwxxoQF+/uISAXi69lWpGGEZ6y1ib7XFwDGmDbA5UBbYDDwsjHGa4zxAi8B5wBtgCt85wL8y1dWc2A/cH2wv4iIVCAhvrAlK6DFlrVhhKHAe9baDGvtemAt0NX3WmutXWetzQTeA4YaYwxwJvC+7/qJwIWl0G4RqSi8zi/HFalne4sx5hdjzJvGmBq+Y7HA5nznbPEdK+p4LeCAtTa7wHERkVMTUs6GEYwxM40xywt5DQXGAc2ARGA78N+8ywopyp7C8aLaNMYYs9AYs3D37t0n9X1EpJLwzbMNsYEdRggJaGn5WGvP8uc8Y8xrwGe+H7cAjfJ9HAds8/25sON7gOrGmBBf7zb/+YW1aTwwHiApKanIUBaRSswXtuWmZ1scY0yDfD9eBCz3/XkacLkxJtwYkwA0B34GFgDNfTMPwnBuok2z1lrgO+AS3/UjgU+C8R1EpILyuvMEmWs92xP4tzEmEedX/g3AWABr7a/GmKnACiAbuNla5wFlY8wtwHTAC7xprf3VV9Z9wHvGmH8AS4A3gvlFRKSCCXFukIUEeDZCqYSttfbqYj57HHi8kONfAF8UcnwdzmwFEZGS8+aN2Waf4MSTU9amfomIlK68nm2Ab5ApbEVE8vPm3SBT2IqIuCdv6lcFf4JMRKR0ebxgvOrZioi4LiRcY7YiIq7zhhKCZiOIiLjLG4ZXU79ERFzmDdMwgoiI6zSMICISBN4wPUEmIuI6T6iGEUREXOcNxUtOQItU2IqIFOQN1WwEERHXeULx6gaZiIjLvCG6QSYi4jqPhhFERNynebYiIkHgCcFrNRtBRMRdXt0gExFxnydUPVsREdd5QvRQg4iI67whgIgHDQAACJdJREFUmo0gIuI6jx7XFRFxn2YjiIgEgScEj3q2IiIu86pnKyLiPvVsRUSCwBNCCDlgbeCKDFhJIiIVhSfEeSM3cEUGrCQRkYrC4wUI6PQvha2ISEHGCVv1bEVE3JQ3jBDAGQkKWxGRgo4OI6hnKyLinqM3yNSzFRFxj69n67Hq2YqIuMdoNoKIiPs8mo0gIuK+vKlfmo0gIuIi9WxFRIJAYSsiEgR5N8g0jCAi4iL1bEVEgkBrI4iIBIGvZ2sUtiIiLjJONOoJMhERN+WFrXq2IiIu0qpfIiJBoBtkIiJBoFW/RESCwDdmq9kIIiJuMnlTv7SVuYiIezQbQUQkCDwKWxER9/l6tt3iqwesSIWtiEhBvjHbM1vUCliRClsRkYJ8PVusbpCJiLjnaNhWgDFbY8ytxphVxphfjTH/znf8AWPMWt9ng/IdH+w7ttYYc3++4wnGmPnGmDXGmCnG/H979xoi11nHcfz7Y+PGooE2jdo0TWxSNi9SKmnclGBtRRFtInVb6YuI0LyoBkoLihQaCUgUX1ipFoQqRLIYtbiGprVBrLcS2r6x6UbTXAxJ115oLjSVegtqkpq/L84TO07PzF6cec6c2d8HDnP2mWfO/v88s/89lznPaDB3LmbWZ/ql2Er6MDACvC8irgbuT+0rgPXA1cBNwHckDUgaAB4E1gIrgE+nvgD3AQ9ExBDwZ+COrMmYWf9Jd5DVvtgCdwJfj4gzABFxKrWPAGMRcSYiXgQmgOvSMhERL0TEWWAMGJEk4CPAw+n124FbMuZhZv1IKh77oNguB25Ih/9PSlqd2hcBrzT0O5baWrVfCvwlIt5oajczm7kunEaY07EtNZH0G+Cykqc2p997CbAGWA3skLQMUEn/oPyfQrTp3yqmjcBGgCVLlrQL38xmszlvh0XDcNH8zm2yY1tqEhEfbfWcpDuBRyIigD2SzgMLKPZMFzd0vQI4kdbL2v8EXCxpTtq7bexfFtNWYCvA8PBw5z7TYWb9Zd5l8LknOrrJqk4j/JTiXCuSlgODFIVzF7Be0lxJS4EhYA/wLDCUPnkwSHERbVcq1ruB29J2NwCPZc3EzGwKurZnO4lRYFTSQeAssCEVzkOSdgB/AN4A7ooovrhd0t3AL4EBYDQiDqVt3QuMSfoa8HtgW95UzMwmp+jgHRJ1Mjw8HOPj41WHYWZ9RtLeiBhubvcdZGZmGbjYmpll4GJrZpaBi62ZWQYutmZmGbjYmpll4GJrZpaBi62ZWQYutmZmGbjYmpll4GJrZpaBi62ZWQYutmZmGbjYmpll4GJrZpbBrJ3PVtJrwMtT6LqA4lsk+oFz6V39lM9sz+W9EfGu5sZZW2ynStJ42UTAdeRcelc/5eNcyvk0gplZBi62ZmYZuNhObmvVAXSQc+ld/ZSPcynhc7ZmZhl4z9bMLAMX2xYk3STpiKQJSZuqjmcmJL0k6YCkfZLGU9t8Sb+W9Hx6vKTqOMtIGpV0StLBhrbS2FX4dhqr/ZJWVRf5W7XIZYuk42ls9kla1/Dcl1IuRyR9vJqoy0laLGm3pMOSDkn6fGqv3di0yaU7YxMRXpoWYAD4I7AMGASeA1ZUHdcM8ngJWNDU9g1gU1rfBNxXdZwtYr8RWAUcnCx2YB3wOCBgDfBM1fFPIZctwD0lfVek99tcYGl6Hw5UnUNDfAuBVWl9HnA0xVy7sWmTS1fGxnu25a4DJiLihYg4C4wBIxXH1CkjwPa0vh24pcJYWoqIp4DXm5pbxT4C/CAKvwUulrQwT6STa5FLKyPAWESciYgXgQmK92NPiIiTEfG7tP534DCwiBqOTZtcWvm/xsbFttwi4JWGn4/RfhB6VQC/krRX0sbU9p6IOAnFmw14d2XRTV+r2Os6XnenQ+vRhtM5tclF0pXAtcAz1HxsmnKBLoyNi205lbTV8WMb10fEKmAtcJekG6sOqEvqOF7fBa4CVgIngW+m9lrkIumdwE7gCxHxt3ZdS9p6Kp+SXLoyNi625Y4Bixt+vgI4UVEsMxYRJ9LjKeBRikOeVy8cxqXHU9VFOG2tYq/deEXEqxHx74g4D3yPNw9Hez4XSW+jKE4PRcQjqbmWY1OWS7fGxsW23LPAkKSlkgaB9cCuimOaFknvkDTvwjrwMeAgRR4bUrcNwGPVRDgjrWLfBdyernyvAf564ZC2VzWdt7yVYmygyGW9pLmSlgJDwJ7c8bUiScA24HBEfKvhqdqNTatcujY2VV8R7NWF4irqUYorjpurjmcG8S+juHL6HHDoQg7ApcATwPPpcX7VsbaI/8cUh3DnKPYo7mgVO8Xh3YNprA4Aw1XHP4Vcfphi3Z/+iBc29N+ccjkCrK06/qZcPkhx6Lwf2JeWdXUcmza5dGVsfAeZmVkGPo1gZpaBi62ZWQYutmZmGbjYmpll4GJrZpaBi62ZWQYutjbrpCn07plG3+OSvjrN3/GQpNcl3TazKK3fuNiaTe6BiPjydF4QEZ+hZncdWnfNqToAsxwkbQZup5i16TVgr6TLgZ83dLsGWBYRL7fZzhaKuUwXAsuBL1LM07oWOA7cHBHnupGD1Zv3bK3vSXo/xfwW1wKfAlZDMVFPRKyMiJUUE47sbFdoG1wFfIJiftMfAbsj4hrgn6nd7C28Z2uzwQ3AoxHxDwBJ/3N4L+l64LOp31Q8HhHnJB2g+FaPX6T2A8CVHYnY+o6Lrc0WpZOApBmetgGfjIjTU9zWGYCIOC/pXLw5wch5/DdlLfg0gs0GTwG3SrooTTt5M/x3LtMdwL0RcbTKAK3/udha34vie6Z+QjGF3k7g6fTUByjO336l4ZtUL68oTOtznmLRrI306YPTEXH/DF77feBnEfFwp+Oy+vGerVl7p4GNM7mpAfgQ8K+uRGW14z1bM7MMvGdrZpaBi62ZWQYutmZmGbjYmpll4GJrZpbBfwAvNhfhljobbwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,8))\n",
    "ax.plot(dz_new, -eh_new, label='original')\n",
    "ax.plot(ds.dz, -ds.z, label = 'initial guess using func_tanh')\n",
    "ax.legend()\n",
    "ax.set_xlabel('dz [m]')\n",
    "ax.set_ylabel('depth [m]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, the inital guesses get pretty close to the original profile. We can do a much better fit by using curve_fit to come up with coefficients for function func_tanh using our initia guesses as input."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimize the parameters to get a better fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "popt, pcov = curve_fit(func_tanh, eh_new, dz_new,p0=[a, b, c ,d])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = 108.77989339824059, b = 140.2328276377167, c = -0.14268724607567607, and d = 984.7613639700135\n"
     ]
    }
   ],
   "source": [
    "print('a = {}, b = {}, c = {}, and d = {}'.format(popt[0],popt[1],popt[2],popt[3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2b4b2b1ece80>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAHSCAYAAABowJlSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3xN9//A8dfn3kyCIImRhBgxEwmxYkSMoDYtqtQqihZdqq1WW/otVe2v06itpbZaVTsooYJEbWolZowgkX3P748baVQSGfdmeT8fjzwq557z+XxOom/nfMb7ozRNQwghRObo8roBQghRkEjQFEKILJCgKYQQWSBBUwghskCCphBCZIEETSGEyAKLvG5ATjk4OGhubm553QwhRCFz6NChW5qmOf73eIEPmm5ubgQHB+d1M4QQhYxS6lJax+X1XAghskCCphBCZIEETSGEyIIC36cp8kZCQgLh4eHExsbmdVOEyBEbGxtcXFywtLTM1PkSNEW2hIeHU6xYMdzc3FBK5XVzhMgWTdO4ffs24eHhVKpUKVPXyOu5yJbY2FhKly4tAVMUaEopSpcunaU3JgmaItskYIrCIKt/jyVoikKvQ4cOREZGZnjOhAkT2LZtW7bKDwwMpFOnTtm6VhQ80qcpCi1N09A0jd9///2p506cODEXWiQKA3nSFAXa119/jYeHBx4eHnzzzTdcvHiRmjVrMnLkSOrVq0dYWBhubm7cunULgEmTJlGjRg0CAgLo06cP06ZNA2DgwIGsXLkSMK4y+/jjj6lXrx6enp6cOnUKgL/++osmTZpQt25dmjRpwunTp/PmpkWekidNkWOfrj/Oiav3TVpmrfLF+bhz7QzPOXToEPPnz+fAgQNomkajRo1o0aIFp0+fZv78+UyfPv2x84ODg1m1ahVHjhwhMTGRevXq4ePjk2bZDg4OHD58mOnTpzNt2jTmzJlDjRo12L17NxYWFmzbto0PPviAVatWmeyeRcEgQVMUWH/++Sfdu3enaNGiAPTo0YM9e/ZQsWJFGjdunOb5Xbt2xdbWFoDOnTunW3aPHj0A8PHxYfXq1QDcu3ePAQMGcPbsWZRSJCQkmPqWRAEgQVPk2NOeCM0lvU0BHwXRzJ6fFmtrawD0ej2JiYkAfPTRR7Rs2ZI1a9Zw8eJF/P39s9ZgUShIn6YosPz8/Pjtt994+PAh0dHRrFmzhubNm6d7frNmzVi/fj2xsbFERUWxcePGLNV37949nJ2dAViwYEFOmi4KsHwXNJVS7ZVSp5VS55RS7+V1e0T+Va9ePQYOHEjDhg1p1KgRQ4YMoWTJkume36BBA7p06YKXlxc9evSgfv36lChRItP1vfvuu7z//vs0bdqUpKQkU9yCKIBUftr3XCmlB84AAUA4cBDoo2naifSuqV+/vib5NHPfyZMnqVmzZl43I8uioqKws7Pj4cOH+Pn58dNPP1GvXr28bpbIY2n9fVZKHdI0rf5/z81vT5oNgXOapp3XNC0eWAp0NVXhQWFB+M7xxfVrV8ZtG5fxyQYD3LlgqqpFPjFs2DC8vb2pV68ezz//vARMkWX5bSDIGQhL9X040MgUBQeFBdFsXjMMGECDqXunsv/iOT5qMh07GwvsrI1f9kUssS9iBYGfw4GfoO8KqGCSJoh8YMmSJXndBFHA5begmdYi0Cf6D5RSw4BhABUqVMhUwYEXA40B81EtGuwOX03vxeUolvTcY+cuGtwQv3oD4Nhq+LkbvLgYqrTK0o0IIQqn/PZ6Hg64pvreBbj635M0TftJ07T6mqbVd3R8Yt+jNPm7+aNSx2Rl/LprNZ3eLU4z62Uf3mxTDYAb92PB3hUG/wGlqsCS3nByffbvSghRaOS3oHkQcFdKVVJKWQEvAutMUbCvqy8zO818PHACGhrT/hpL8eIX6FHP+fGL7Jxg4Hoo5w3LB0DIr6ZoihCiAMtXQVPTtETgdWAzcBJYrmnacVOVP8xnGHsH76WWQ63HjhswMGTdEA5dO/DkRbYl4eU1UKk5/Dbc2M8phHhm5augCaBp2u+aplXTNK2Kpmn/M3X5vq6+zOkyB516/NZP3DrBC6sCeKDf9ORF1nbQZxnU6ASbxsLuLyEfTdUSQuSefBc0c4Ovqy8zOs544lU9SUvijuV0tl1e/ORFljbQcyHUeRF2fAZbPjROSxJ5pkmTJk89Z8iQIZw4YZzm+/nnn2f5ejs7u+w1Lp/JTE7RzBo7diy1a9dm7NixJikvPSEhIZlK65cRc/z+8tXk9uzIyeT2nw79xMiNI0nSUq3u0IyZnGd2mskwn2FPXmQwwB/j4K+foHYP6DbDGFCfMQVxcrudnR1RUVFmv6awK168OBERESnr881lwYIFBAcH88MPP2S7jMz+/gry5PZcNcxnGHsG7Xm8j1MZB4eGbxjOT4fS6L/U6eC5qRAwEY6vhl96wMM7udfo/GjTezC/o2m/Nj19Be2jp4jAwED8/f154YUXqFGjBn379k1JzuHv709wcDDvvfceMTExeHt707dv38euj4qKonXr1in5M9euXfvUutPLy/moPoBbt27h5uYGQFJSEmPHjqVBgwbUqVOHWbNmAXDt2jX8/Pzw9vbGw8ODPXv2kJSUxMCBA/Hw8MDT05P/+7//e6L+1Pk/U99LWuUBKTlFH+UbHTp0KLVr16Zt27bExMQAcPDgQerUqYOvry9jx47Fw8PjiXq7dOlCdHQ0jRo1YtmyZem2I6PfycGDB2nSpAleXl40bNiQBw8ePFFPfHw8EyZMYNmyZXh7e7Ns2bJ085kuWLCAHj160L59e9zd3Xn33XcfK2v8+PF4eXnRuHFjbty4keHvNTOe6aAJ//ZxWuqSt+/UHv0ng8CpFDQdA8/PhfCDMK8d3L2Ue40WTzhy5AjffPMNJ06c4Pz58+zdu/exz6dMmYKtrS0hISEsXvx494uNjQ1r1qzh8OHD7Ny5k7fffjvDjEip83KuXr2azLzpzJ07lxIlSnDw4EEOHjzI7NmzuXDhAkuWLKFdu3aEhIQQGhqKt7c3ISEhXLlyhWPHjvH3338zaNCgTP8c0irvv86ePctrr73G8ePHsbe3T8kJOmjQIGbOnElQUBB6vT7N8tetW5fyc+zdu3eGbUnrdxIfH0/v3r359ttvCQ0NZdu2bSmp+lKzsrJi4sSJ9O7dO6WuR/lMjxw5wsSJE/nggw9Szg8JCWHZsmX8/fffLFu2jLAw4xqZ6OhoGjduTGhoKH5+fsyePTvTP8v05LfJ7XnC19WXXQN30X/NYM7dOZVy/FHg3HRuE+82eRdfV9/HL/R8AYqVhaUvwZw20Hc5lK+by63PB56bktctoGHDhri4uADg7e3NxYsXadasWaau1TSNDz74gN27d6PT6bhy5Qo3btygbNmyaZ6flbycj2zZsoWjR4+mPJXdu3ePs2fP0qBBAwYPHkxCQgLdunXD29ubypUrc/78eUaNGkXHjh1p27Ztpu4DSLO8/6pUqVLKcR8fHy5evEhkZCQPHjxI6ed96aWX2LBhQ6brTUtav5MSJUpQrlw5GjRoABhf9TMro3ymrVu3Tkm+UqtWLS5duoSrqytWVlYp+zf5+PiwdevWHN0TyJNmCl9XX+Z2mQPq8X9hNTR+O/UbLRa0ICgs6MkL3ZrBK1vBwgbmd4AzW3KpxSK11P1rqXNgZsbixYuJiIjg0KFDhISEUKZMmQy3dM3oKdTCwgJD8gBh6jI0TeP7778nJCSEkJAQLly4QNu2bfHz82P37t04Ozvz8ssvs2jRIkqWLEloaCj+/v78+OOPDBkyJMN6NE0jPj4eIM3y/iutn1V2xzbSa0dG9WR3F9NH+UyPHTuWkuIvo7oALC0tU+rL6t+L9EjQTMXPrSnDaizCNrHxEyPrCYYEpu6dmvaFjtVhyFZwcIdfX4TgebnQWpFVlpaWaWZbv3fvHk5OTlhaWrJz504uXcq4qyWjvJxubm4cOnQI4LG+vnbt2jFjxoyU+s+cOUN0dDSXLl3CycmJoUOH8sorr3D48GFu3bqFwWDg+eefZ9KkSRw+fPiJNqSuZ+3atSnlplVeZpQsWZJixYqxf/9+AJYuXZqp69JrR3pq1KjB1atXOXjwIAAPHjxIN5AVK1bssf7O/JLPVILmf/xf9160cppK6cTXngicv53+je7Luqf9xFmsLAz83bhGfcOb8Mf7kJTzf9WE6QwbNow6deqkDAQ90rdvX4KDg6lfvz6LFy+mRo0aGZaTUV7Od955hxkzZtCkSZOUzdzAOPWpVq1a1KtXDw8PD1599VUSExMJDAzE29ubunXrsmrVKsaMGcOVK1fw9/fH29ubgQMHMnny5CfaMHToUHbt2kXDhg05cOBASrb6tMrLrLlz5zJs2DB8fX3RNC1TuUbTa0d6rKysWLZsGaNGjcLLy4uAgIB0n+pbtmzJiRMnUgaC8ks+02d6ylF67j1M4KU5+zkUsYKbFj+g/SdniKXOklfqvkJ/r/5P9nMmJcKW8XBgJlRpDS/MA1t7k7YvPyiIU45MqTDm5Xx0T2AcOLt27RrffvttHrcqd8iUoxwqUcSSRYMbUrNED8okjkL958eUYEhg5qGZ+C3we3J0XW8Bz30Bnb+DC7thTmu4dTYXWy9yQ2HMy7lx48bHpip9+OGHed2kfEmeNDNw80EsrywIZt/1pdy1monGk68EOqVjRscZaU+Ev7QPlvUzPn32nAdV25ilnXnhWX/SFKazefNmxo17PCl4pUqVWLNmTa61IStPmhI0nyI2IYmP1x5n0eEtWBf/k2rlk9hxcdNjq4j0Ss/QekPTfl2/e8k4JenmCWj7P2g8wjjPs4CToCkKE3k9NyEbSz1fvFCHOS/2pUTcSC6dGcmL7p8+lvAjSUti1qFZtF7U+slBopIVYfBmqN4BNr8P616HxLhcvgshhKlI0MykjnXKsfVNP1pUc+TPUG/q2r2Dhe7ftQEaGrGJsXwS+MmTgdPaDnr9DH7vwpFfYEEnuP9EbmUhRAEgQTMLnIrbMOtlH77u5UVMZEsqJE2lhfNLWOmsAGPg3HJ+C/4L/RmxYcTjwVOng1bjoecCuHEcZvnBhT15cyNCiGyToJlFSil61HNh4+jmNHFtwsVzL1HL8ivqOvmlnBOfFJ/+63rt7jB0B9jYw6KusPc7yc0pRAEiQTOb3ByK8vMrjVgwqAEOVp5cC+uMXlmnTIjX0IhJjGHCzglPBk6nGsbAWaMjbP0IlveH2Pt5cBcF18WLF9PMwpNVgYGB7Nu3L83P4uLiaNOmTcrk6oxyc2ZXcHAwo0ePNklZ+UFERASNGjWibt267NmzJyWPZ2RkJNOnT8/r5pmEBM0c8q/uxO9jmjOhbXecEz7H3vAcFsoq5fNtF7bRcmHLJwOnTXHotQjafganNsLsVnDzZC63PncFhQUxec/ktFdU5ZGMguaRI0dISEhIybIzZ84catUyphE0VdCsX78+3333nUnKyglTrMkG2L59OzVq1ODIkSM0b96c33//HXt7+0IVNGXKkQldvv2QzzaeYP2pQB7aLOW+dphHuea8ynjRo2YPAioHPDkt6eJeWDEQ4qOh6/fg8Xyutz2rUk/ReOOPNwi5HpLh+ffi7nH0xlEMmgGd0lGnTB1KWKe/TM+7rDfftP8m3c8vXrxI+/btadSoEUeOHKFatWosWrSIIkWKcOjQId566y2ioqJwcHBgwYIFlCtXju+++46ZM2diYWFBrVq1mDJlCo0bN0av1+Po6Mj3339P8+bNAbh58yZNmjQhIiKCSpUqsWrVKl555RWmTZvGypUr+fLLL/H09KR27dpPpJpLnfh25cqVbNiwgQULFrBixQo+/fRT9Ho9JUqUYPfu3QQGBjJt2jQ2bNjAJ598wuXLlzl//jyXL1/mjTfeSHkKnTRpEosXL8bV1RUHBwd8fHx45513Hqv3xo0bDB8+nPPnzwMwY8YMypcvT6dOnTh27BgA06ZNIyoqik8++QR/f3+aNGnC3r17adWqFfPnz+f8+fPodDoePnxI9erVU9ry2muvERERQZEiRZg9e3aaS01DQkLo0qULMTExODs7ExQURM2aNQkODub1119n7dq1VK9enYCAAL788ssM/77ktqxMOZLUcCZUoXQRfupfn0OXqvDu2uLsvDMKTSWigNAboYTeCOXzPZ+zc8DOxwOnW1N4dbcxcK4cDGEHjUmOLazSq6rAuRd7D4NmzIZj0Azci72XYdDMjNOnTzN37lyaNm3K4MGDmT59OmPGjGHUqFGsXbsWR0dHli1bxvjx45k3bx5TpkzhwoULWFtbExkZib29PcOHD8fOzu6JAOTk5MScOXNSAlpqU6ZM4YcffiAkJON/KP5r4sSJbN68GWdn53S3njh16hQ7d+7kwYMHVK9enREjRhAaGpqSvzMxMZF69erh4+PzxLWjR4+mRYsWrFmzhqSkJKKiorh7926GbYqMjGTXrl0AHD58mF27dtGyZUvWr19Pu3btsLS0ZNiwYcycORN3d3cOHDjAyJEj2bFjxxNleXt7M3HixDSzrU+ZMoVjx45l+WeWH0nQNAOfiiXZNmooP+x15vs/13ItOpwoiz8AjbikOAauHUivWr3o4N7h3+BZvBwM3ABbJ8D+6RC237huvVTlPL2XzMjoifCRoLAgWi9qTXxSPFZ6Kxb3WPzkE3cWubq60rRpUwD69evHd999R/v27Tl27BgBAQGAMWN6uXLlAFKSdXTr1o1u3brlqO7saNq0KQMHDqRXr1706NEjzXM6duyItbU11tbWODk5cePGjUzn79yxY0dKKrhHT7NPC5qpEwn37t2bZcuW0bJlS5YuXcrIkSOJiopi37599OzZM+W8uLhne56x9GmaiVKKUc06curdWXzYcjg6rEDTATrO3D7DZ3s+w3+BP/sup+pP01tC+8nQ+xe4cx5mtYBjq/PsHkzJ19WX7f23M6nlJLb3357jgAk8kZdRKYWmadSuXTslb+Xff//Nli3GHKcbN27ktdde49ChQ/j4+JisHy+jdqXO4DNz5kw+++wzwsLC8Pb25vbt209ca8pcl/B4vsv/tgd4LCtRly5d2LRpE3fu3OHQoUO0atUKg8GAvb19ys8zJCSEkycLd9/700jQNDOdTjGuVRcCB26nV7W3cVIdQDP+TxVviKfvmr4sPrr48QGSmp1h+J/GPJ0rB8H6NyAhJg/vwjR8XX15v/n7JgmYAJcvXyYoyPgz+/XXX2nWrBnVq1cnIiIi5XhCQgLHjx/HYDAQFhZGy5YtmTp1KpGRkURFRT2RszGz0svNCVCmTBlOnjyJwWB4bP30P//8Q6NGjZg4cSIODg4pWzI8TUb5O1Nr3bo1M2bMAIxP2Pfv36dMmTLcvHmT27dvExcXl2E2djs7Oxo2bMiYMWPo1KkTer2e4sWLU6lSJVasWAEYEw2HhoZmqt2pZffnnB9J0MwlzSs2ZdlLU1n28jgs9daADjQLwiKv029NP8bvGP/4vE77CjBoEzR9Aw7Nh9mtIeJ0nt5DflOzZk0WLlxInTp1uHPnDiNGjMDKyoqVK1cybtw4vLy88Pb2Zt++fSQlJdGvXz88PT2pW7cub775Jvb29nTu3Jk1a9bg7e2dsglZZqSXmxOM/XedOnWiVatWKV0DYNz61tPTEw8PD/z8/PDy8spUXRnl70zt22+/ZefOnXh6euLj48Px48extLRkwoQJNGrUiE6dOj01V2jv3r355ZdfHnttX7x4MXPnzsXLy4vatWtnauO5/ypdujRNmzbFw8PD7Fv/mpuMnueBoLAgNp/bTuRdd5Yc/Y0IlvIo33GvWr0Y02gMuy7twt/N3/hUdnYbrHkVEh5Ch2ng/VKeJ/2QhB25qzDm78xPZPQ8n/N19U15Re3kVZ4OS34jwRAPmsbyE8tZcWIlSoG13trY/+fexvi6vnoorB0JF3ZBx6/Aulge34nILcOGDePEiRPExsYyYMAACZh5SIJmHmtTpTm7Bu5gyz87iImqxqzD3xJp2IsGxCTGsv7MemOALV4O+q+F3dNg1xTj1sE95oDLk1NPROGzZMmSvG7CY/73v/+l9HM+0rNnT8aPH59HLco98nqez+y6+CcBP7chwRAHGuiUNS/Xeov+Pm05cCXI+MqeZDC+rt+/Ci3fh2ZvgS7tfarNRV7PRWEir+cFWAu3ZuwauJPtF3ZyO9KeX44uZuGJySw8MRmFwsbCxvjKPvxP2Pg27PgMzm2H7rOMuTtzUU62YxUiv8jqg6M8aeZzSQaN5xa9xNZL/26p2rz8i/ze/2fsrC3g6HLY8JZxYKjj11CnZwalmc6FCxcoVqwYpUuXlsApCixN07h9+zYPHjygUqVKj30m210UYI9W08QmxqbsjFlMa8TQOh8yLqANTonXYfWrxlVEnj2Ng0Q2OVui+DQJCQmEh4enu/2qEAWFjY0NLi4uWFpaPnZcgmYBFxQWRODFQBq7NGb18Z38dORr4pNiKG4IoE2lDlQqe4vnE27je+RXKO4MPWZBxSZ53WwhCiwJmoXMrYe3GLv5YxYenYWWvMmbXmfNwkZf8tKpuajIS9DsTWjxXqFK/CFEbpGN1QoZhyIOzO/+I2/7vmmcGK8gSYvjlb1zeUF9yuUK3WDPV8Y8ndeP5XVzhSg0JGgWcD1q9sDWwha90qNTOuL0oay/OwDvf/S8ZTmWh3evov3kD3/+Hxie3LddCJE18npeCDzq7/R380ev0zN+x4dsO78VG50DxeIa08HiJK9yDZ9yvli9MAtKV8nrJguR70mf5jMm8GIgozeN5u+bfwOg03Rs0OxppbPgdpMPKd/m9Txfvy5EfiZ9ms8Yfzd/XvR4EV3yr9igDDxvkcDHhuI47h1P6JRW7A4OxWAo2P9oCpHbJGgWYi3dWmJtYY1e6bHSW1GxlDNf6M5RwdKCfXHB1F7fjslffMqifRd4GG+ehLxCFDYSNAux1NnSAwcEcmLkCTb02YCLUxVGqijqWdyjTNwX3NrUA8/P+zFi+a9cu1fwkx0LYU7Sp/kM0jSNTec28Ungxxy8GvwolSd6zYIyCZN53rM1rzSrRB0X+zxtpxB5Sfo0RQqlFB3cO3BgyF8M8BqAhnGj4UQScS/yJYdPHqXLD3vpNSuILcevkyT9nkKkkKD5DFNK8arPq9ha2KJTOnRKEZh4k1O6wbSpMJNLt28y7OdDtP4qkEVBF6XfUwjk9Vzw+DzPpLsX+fz3UWyKu01xnSUd3fpzN7IbJ68oStha8lKjCgzwdaNsCZu8brYQZiXzNEXmGZI4suUDPj/wLauIw1ZvTddqA9E/7MLe0xo6pejsVZ5XmlXCw9m82ZSEyCsSNEXW3TrLyZUD+OJ6ML+oRHQ6C/wrtCXqoQPXrnmhxVejceVSDGlWmVY1nNDpZLK8KDwkaIrsMSTBgZlc3PYJb2vRrNaMU5L0Sk/f6v/j7Pn6XL0XS2WHogxuVonn67lga5W7W28IYQ4yei6yR6cH39dwG7mf+iUqoE/+NzZJS2LRqfewLT+Zga3uUdRaz4e/HcN3ynambT7NzQeSnFgUThI0ReaUroJ/97lY6S3RA7bAaxVacOr2ST4N6ku41WheaXuV+hWL82PgOZpN2ck7K0I5df1+XrdcCJOS13ORJUFhQQSeWov/pQP4hh8mzqUhS6q3ZurRhZy6dQo3ezcGer5OTGRz1h65Q0xCEs3dHXilWSVaVHOU/YREgSF9msK0NA2OLoNN4yAhBoP/ODY4VOaLoGnsC9tHadvSDKk7gpJaZ1YdvM/NB3G4O9kxpHkluno7Y2Mp/Z4if5OgKczjwQ34/R04uQ7KeUPXH9gbf58v9n7B+jPrsbWwpX3VDiTFO3DjphfXIyrgYGfFy43d6Ne4AqXtrPP6DoRIkwRNYV4n1hr3YY+5C03fAL+xnIg8z9itY/n97O+AccT93YZfc+1aY3aejsDaQkev+q4MbV6ZCqWL5PENCPE4GT0X5lWrK7z2F9TpDXumwazm1Iq5TzPXZuiV8VU8SUti8oExXFDv8kGPaLp6lWPpwcv4T9vJ60sO83f4vTy+CSGeToKmMJ0ipaDbdOi3GhJjYV57/K8dw0pvhV7psbWw5fUGr/PP3X94dVNv1t/sw4gOlxjYzJldpyPo/MOf9J2znz1nIyjob0Ci8JLXc2EecVGw4zM4MJOgoiUJdG+Jv89QfF19SUhKYPnx5UwLmkbI9RDKFC3D0HojKWHowPIDkdx8EEft8sV5tUUVOniUxUIv/7aL3Cd9miJvhP0F60ZBxCnjq3u7yVC0NGDM67njwg6mBU3jj3N/UMSyCAO8BlG9aG/WHkrifEQ0rqVsGdq8Mj19XGWlkchVZunTVEr1VEodV0oZlFL1//PZ+0qpc0qp00qpdqmOt08+dk4p9V6q45WUUgeUUmeVUsuUUlY5aZvIJ1wbwqu7ocV7cGw1/NgQ/l4JmoZSitaVW7Op7yaODj9Kr9q9mHP4J97a7U8p1x95s4PCwc6aCWuP02TKdr7Zdoa70fF5fUfiGZejJ02lVE3AAMwC3tE0LTj5eC3gV6AhUB7YBlRLvuwMEACEAweBPpqmnVBKLQdWa5q2VCk1EwjVNG3G09ogT5oFyI0TsO51uHIIqrWHjl9DCefHTrn64CrfH/ieGcEzuBd3j2auzfBybMX+fyK4crMy9noPejdw5ZVmlXAtJSPuwnzM+nqulArk8aD5PoCmaZOTv98MfJJ8+ieaprVLfR4wBYgAymqalqiU8k19XkYkaBYwyQlA2PEZKD0EfAI+g0H3+EvPg7gHzDsyjyl/TuF69HUALHSWdHOZzpGz5dGArl7lGeFfBfcyxXL/PkShl9tTjpyBsFTfhycfS+94aSBS07TE/xwXhU1yAhBGBoGLj3Fu54KOcOvsY6cVsy7GmMZjeL3h6ynbECcaEth07Q26twiiZ8MSbDp2nYD/282wRcGEhEXmxd2IZ9BTg6ZSaptS6lgaX10zuiyNY1o2jqfXpmFKqWClVHBERETGNyDyp5Ju8PJv0HU63DwOM5rC7mmQlPDYaa0qtUrZhthab41XGS+m7f+M744H4Ouzhj5NLNl//jbdftxL3zn72XvulkxXEpqXQGIAACAASURBVGZl8bQTNE1rk41ywwHXVN+7AFeT/5zW8VuAvVLKIvlpM/X5abXpJ+AnML6eZ6N9Ij9QCur2haptYNNY2DEJjv8GXb4D53rAv9sQP9qOw9fVlxMRJ/g66Gt+PrqAhKTZdKrelZp2L7HjqDV95xzAy9Wekf5VCKhZRhIjC5MzV59mbWAJ/w4EbQfcMT5RngFaA1cwDgS9pGnacaXUCmBVqoGgo5qmTX9a3dKnWYic3GB8XY++aXyF9/8ArNIf7LkedT1l0Ohu7F0au/jSyHEAwScrE343HncnO0a2rELnOuVlrqfIMrMMBCmlugPfA45AJBCSapBnPDAYSATe0DRtU/LxDsA3gB6Yp2na/5KPVwaWAqWAI0A/TdPintYGCZqFTEwkbJ0AhxdCyUrGp85KfhleEhUfxfwj8/l6/9dcjLyIeyl3WrsM5uwFH87djMelpC2v+lWmZ31Xya4kMk0mt4uC5cJuWDca7l6AegMgYCLY2md4SaIhkTUn1/Dlvi85ePUgDkUceM5tALdu+HMi3Djnc0jzSvRtVIFiNpa5dCOioJKgKQqe+IcQOBmCfgC7MsZ5nTU6PPUyTdPYc3kP0/ZNY/2Z9dhY2NCwbCsiIktw93Y9HK086e/rxqCmbpKaTqRLgqYouK4chrWvG0fZa/eA56aCnWOmLj116xTvbn2X9WfWA6BTOpqXHs+l8EbYWuoZ1NSNYc2rUKKIPHmKx0lqOFFwOdeDYYHQ8kM4tQF+bAChy4zZ45+ihkMNfF18U9LTGTQDu25NonzVr6hW8SI/7DxHs6k7+H77WaLiEp9SmhASNEVBYWEFLcbCq3ugtDusGQaLe0Jk2FMv9Xfzfyw93cj6IzkfeZK14SMoUeFjypQJZdrWU/hN3cns3eeJTUjKhRsSBZW8nouCx5AEf82G7RONcz3bfAL1X3liKWZqQWFBj831jE2MZVHoIqbunco/d//BrYQ7ZXW9uHq1PmWK2TGqVVV6NXDF2kJG259V0qcpCp+7l2DDG/DPDqjYFLp8D6WrZKmIJEMSq06uYsqfUzhy/QiORcpRXteTOxHNcbUvxZjW7vSo5yzzPJ9BEjRF4aRpELIY/vgAkuKh9UfQaLhxjXuWitHY8s8WpuydQuDFQIpZ2VNO342HdwKo6lCeN9q407lOeVlh9AyRoCkKt/vXYMObcGYTuDSArj+CY/VsFbU/fD9f7P2C3079hrXeljL6jiTd70Btpyq81bYabWuVkf3bnwESNEXhp2nGBMebxkJ8NPi/B03GgP6pKRbSdCLiBFP3TmXx34vRNHDUt0ZFdaW+sydvBVSjRTVHCZ6FmARN8eyIumnci/3E2uS92H+Esh7ZLu7yvct8HfQ1sw/P5mHCQ0rpmkJMPcqXjmd0s24MbdTehI0X+YUETfHsOf6bMXjG3IXm70Dzt41Tl7Lp1sNb/PDXD3wd9DUP4h8YD2qWdCg7g2979KKqkyRDLkxkcrt49tTuZtyLvXYP2DUFZreEa0ezXZxDEQc+8f+Et3zfQj1KAasS2Hb9Y1p8O5cJa49xR/YwKvQkaIrCrUgpeH429FkK0RHGwBn4xRPJjrOiXZV22FjYoFd6LHWWWFs/4KrVm3x16BUafTmbOXvOE59oMOFNiPxEXs/Fs+PhHdg0Dv5eDmXrQPeZUKZ2topKPVnew8mD7//6nql7v+ReXCS2Sb7UKjqY/3XqJCPtBZj0aQrxyMn1sP4NiL1nHGFv+ka2R9hTuxd7j2/2f8OX+74iOuEBRRKb4Vd2ONO6d6J2+RImaLjITRI0hUgt+pZxkOj4GihfF7rNBKcaJin6bsxdpu37iq+DviE26SFFk/zoXf0NJnduj1NxG5PUIcxPgqYQaTm22rjFRnwUtBwPTUZleTVRem4/vM3/dn/Bjwd/ID4pjhJaS0Y3HMcHbVtJBvkCQIKmEOmJumlcTXRqg3E1UbeZ4FDVZMXfjL7JR9s/Z17ITBINCTjqA5jUagLDmvhKf2c+JlOOhEiPnRP0/gWen2vcf31Wczg4J1P5OjPDqagTs7p8Q9hbF+ldcyi3DTsZvtWPSl90ZcPx7E+BEnlDnjSFSO3+VVj7mjFzUtU20OUHKF7OpFWE37vC0DUfsvnSYjRNw8O+O280e4GbMf+kpK4TeU9ez4XILE0zPmlu+QgsbaDT/0Ht7iav5nTEBQat+oCg6yuAJFBgY2HLjv7bJXDmA/J6LkRmKQUNh8LwPcZthFcMhFVDjdsLm1B1x0rsG/4rI+qP5tECo9jEGMZv+x9JBsken19J0BQiPQ7u8MoW8H8fjq2CGU3gfKDJq3nZqye2Frbo0AGKnZc34jrNk23/mL4ukXMSNIXIiN7SOAF+yFawLAKLusKm9yAhxmRV+Lr6sr3/dj5r9Rlb+wbSrcIX3IyOIOCXlrRZ0IOwe0/fB0nkHunTFCKz4h/Ctk/gr1ngVMs42l6mllmqCjx9mYGrPuRy/DL0Oh1jm7zLRy3ew9bS1iz1iSdJn6YQOWVVBDpMhb6r/k3+8ddsk01NSs2/egVOjZ3Pu3V/xzrRh8l7J+L2TTVWnlhJQX/QKegkaAqRVe5tYMQ+cGtuXIr5ax+Ivm3yamws9Uzp2pr9r26kYdFvuBtlQc8VPWk2z5+jN2R+Z16RoClEdtg5Qd8V0P4L+Ge7cZDon51mqcrDuQR/vjmKr/w34Zj0GgfCj+A9sy4jNo7g9kPTB2uRMQmaQmSXUtB4OAzdATYl4OduxrmdiaZPRGyp1zGqZQ2CRk2ho+Ny7BI7MCt4NlW+c+f7A9+TaEg0eZ0ibRI0hcipsp4wLBDqD4Z938HcNnDrnFmqquJox5rh7fix4w9USvqBhNiKjP5jNN4zvdl+frtZ6hSPk9FzIUzp5AZY9zokxhlXEnm9aLaqrkTG8MHqo2w6t45o2/k8NFyje43uTGs7jRtRN1KSJMvqouyRZZRC5Jb7V2HVELi0F+q+DB2+BDNNFdI0jbUhV/l43RGuJK7kgdVyNBLR0NA0DSu9FdtlWWa2yJQjIXJL8fLQf51x98sjP8Ps1sbsSWaglKJbXWd2vB3AS7XG4Bg9gyK6iiQaEknSkohPiifwYqBZ6n5WSdAUwhz0FtB6gnFOZ9R1mNUCjq4wW3Wl7az5rk9dfnyxDfaxI1FYApCkJXHlwRUMmmz0ZioSNIUwJ/c28OoeKFcHVg+BdaNNugTzv7p6O7P19SH4FPmGEgm9qVq8AT8e/JGAnwMIvx9utnqfJRI0hTC3Es4wYAM0exMOL4Q55htdB3AvU4zA0cPoX3sc8Tcm0KTUeA6EH6DOjDqsPLHSbPU+KyRoCpEb9BbQ5hPou9I4UPRTC+P+RGZS1NqCb1/05rNunkRcb0o1puNsV5meK3oyaO0gHsQ9MFvdhZ0ETSFyk3uAMU+nUy1YOQi2fAhJ5pmYrpTi5cYVWTHcF1udCzFXJ9Cl8ussCl2E9yxvgsKCzFJvYSdBU4jcVsIFBm6EBkNg3/fwS3fjlsJm4uVqz8bRzfBzL0fo8fZ0c5lNksFA8/nN+TTwU1lNlEUSNIXICxZW0PEr6DodLh+An/zh6hGzVWdfxIo5/esztl11jpwtQ8XE7+lY9QU+2fUJfvP9OH/3vNnqLmwkaAqRl+r2hcF/GNPLzW0HRxabrSqdTvFay6r8MqQR0bHWnDk5kLd8fuRExAm8ZnqxMGShpJ3LBAmaQuQ153rw6i6o0AjWjoSNb5sl6ccjTao4sHF0czydS7Dqz4r0q7ScumXrMXDtQHqv7M2dmDtmq7swkKApRH5Q1AH6rYEmo4w7YS7sBA+um626MsVtWDK0Ea+2qMyGIwkUfzCRcb6fsubUGurMqMPOC+ZJc1cYyNpzIfKbY6tg7etgYw99foXy3matbuuJG7y1PAQFDA9QfHdkNGdun6GPZx9qlK5Bm8ptnsm165KwQ4iC5PrfsORFiLkDPX6Cmp3NWt3l2w8ZueQQx67c57VWLuy98zFrT68FwMbChh39dzxzgVMSdghRkJT1NCY3dqoFy/rBnq/NshfRIxVKF2Hl8Cb0qOvMjzvCiYmqlLylMMQmxrLq5Cqz1V3QSNAUIr8qVgYGbgCP52H7p/DbCGOeTjOxsdQzracX/RpX4Oh5Z3TKCp3SAzDvyDxCr4eare6CRIKmEPmZpa1xq2D/9yH0V+O+62acCK/TKSZ19eD1ph1wiJ1EU8eR/NxtMUWtitJiQQv+vPyn2eouKCRoCpHfKQX+78EL84wT4Ge3gpunzFid4oMONRnbsguXL7XjrxM1CRywhzJ2ZQj4OYCNZzaare6CQIKmEAWFx/Mw8HdIjIW5AXA+0GxVKaV4K6Aa77avzrrQq0zdeIvtL++itmNtui7tyi9HfzFb3fmdBE0hChIXH+MAUQkX+OUF+Nu8qd5G+lfl48612Hz8BuNXXeb3l7biV9GPl9e8zHcHvjNr3fmVBE0hCpoSLjBoE7g2hFWvQNB0s1Y3qGklJvfwZPfZCMYsOc2KF9bRrUY3xvwxhgk7JzxzSy8t8roBQohssLWHfqth9VDY/D48uAZtPgWdeZ6D+jSsgI2ljndWHGXYoqPMHbCEUjavM2n3JG49vMX3z32PXqc3S935jQRNIQoqSxvouQA2jTPut/7gOnT90ZhByQy613XBxkLPqF+PMGDuIRYOmo5DEQem7pvKnZg7LOq+CCu9eerOTyRoClGQ6fTGLYKLl4PtEyE6Anr/DNbFzFLdc57l+MlSx/BfDtNn9gF+GTKR0kVKM27bOO7G3mV1r9UUtSpqlrrzC+nTFKKgU8q4XXDX6XBhNyzoCFERZquuVY0yzB/YgMt3HtJ7VhAve4xiTuc5bDu/jTY/tyn0WZIkaApRWNTtC32WQsQZWNAB7l8zW1VNqzqw6JWG3HwQR69ZQbR168PKnis5fO0wfvP9WHdqHZP3TC6UW2rkKGgqpb5USp1SSh1VSq1RStmn+ux9pdQ5pdRppVS7VMfbJx87p5R6L9XxSkqpA0qps0qpZUqpwt85IoSpVWsLL682bt42/zmIvGy2qhq4lWLxkEbcj0mk16wgvBwC+KPvH5y/e55uy7rx4c4Pab2odaELnDl90twKeGiaVgc4A7wPoJSqBbwI1AbaA9OVUnqllB74EXgOqAX0ST4X4Avg/zRNcwfuAq/ksG1CPJsqNoH+a40ZkuZ3gDvm28rCy9WepcMaE59ooP+8v/ByasqguoPQ0DBoBuKT4gm8GGi2+vNCjoKmpmlbNE17tCvTfsAl+c9dgaWapsVpmnYBOAc0TP46p2naeU3T4oGlQFellAJaAY9m6i4EuuWkbUI801zqw4D1EB9tDJwRZ8xWVc1yxZkzoD43H8Qx4pdDvFj7Jaz11gBoaPhV9DNb3XnBlH2ag4FNyX92BsJSfRaefCy946WByFQB+NHxNCmlhimlgpVSwRER5uvwFqJAK+dl3PXSkGTs47xx3GxV1a1Qkik9PDlw4Q6bjxRn54CddK7WGYNmYMeFHWarNy88NWgqpbYppY6l8dU11TnjgUTg0a5QKo2itGwcT5OmaT9pmlZf07T6jo6OT7sFIZ5dZWrBoN9BZ2kcVTfjjpc96rnwql9lftl/mXNXyrP2xbX0q9OPjwM/5o9zf5it3tz21KCpaVobTdM80vhaC6CUGgB0Avpq/66nCgdcUxXjAlzN4PgtwF4pZfGf40KInHJwNwZO62KwsCtcDTFbVe+2r0HL6o58uu44+8/fYVanWXiW8eSlVS9x4e4Fs9Wbm3I6et4eGAd00TTtYaqP1gEvKqWslVKVAHfgL+Ag4J48Um6FcbBoXXKw3Qm8kHz9AGBtTtomhEilVCVjhiSb4vBzd7h50izV6HWKb/vUpWLpIoxcfIjbD2B1r9UYNAPPL3+emIQYs9Sbm3Lap/kDUAzYqpQKUUrNBNA07TiwHDgB/AG8pmlaUnKf5evAZuAksDz5XDAG37eUUucw9nHOzWHbhBCp2bvCgHWgt4KFXeDWObNUU9zGkjkDGpBk0BiyMJgyRSvyS49fOHL9CCN/H1ngE3zIxmpCPGsiThtH1C2sja/tJd3MUs2esxEMmPcXbWqWYWY/Hz7d9QkTd09kZseZvFr/VbPUaUqysZoQwsixOvT/zTgdaWEXuHfFLNU0d3fkw4612HLiBt9sO8OEFhNoX7U9ozaN4kD4AbPUmRskaArxLCrraUwt9/AOLOoCUTfNUs2gpm70qu/CdzvOsenYDRb3WIxLcRdeWPECN6PNU6e5SdAU4lnl4gN9VxiXXC7qCjF3TV6FUopJ3TzwqViSd1aEcvWOnlW9VnHr4S1eXPkiiYbEpxeSz0jQFOJZVtEXXlwCt87C0r5m2SLY2kLPzH4+lCxixbBFwbjY1WJmx5nsvLiT8dvHm7w+c5OgKcSzrkpL6DYDLu2FNa+CwWDyKhyLWTO7f33uPIxn+C+HeNGjHyPqj2DqvqmsOrHK5PWZkwRNIQTU6WncLuP4Gtj6kVmq8HAuwVc9vTl06S4f/XaMr9t+TSPnRry85mXe2vxWgcmGJEFTCGHUdAw0HAZBP8D+GWapomOdcoxuVZXlweGsD43g/WbvE5MYw//t/78Ck0ZOgqYQwkgpaD8FanSCP96H47+ZpZo32lSjfsWSfPHHKQ5f+xtdchiKS4wrEGnkJGgKIf6l08Pzc4zbA68eBpf2mb4KneKTLrW5HR3P9ZuVsbb4N42cr4uvyeszNQmaQojHWdoat82wd4Vl/cyS/d3DuQS967uyPdSehV3WM6TeEDQ0gq/l/9V9EjSFEE8qUgr6LIOkRFj6knH1kIm93bY6tpZ6fg8uxuzOs3mu6nN8vudzImMjTV6XKUnQFEKkzaEqvDDPmLx47Wtg4jwVjsWsGd3anZ2nI9h56iaTW08mMjaSL/78wqT1mJoETSFE+tzbQOuPjVOR/vza5MUPaOJGZYeiTNpwgpoOnrzk+RLfHviWK/fNsx7eFCRoCiEy1nQMeLwA2yfBmc0mLdrKQsdHnWpx/lY0i4IuMqnlJBINiXy661OT1mNKEjSFEBlTCrp8D+XqwKohJt+krWUNJ1pWd+TbbWexsyjPiPojmHdkHqdvnTZpPaYiQVMI8XRWRaD3YmMC46UvQVyUSYv/sFMtYhKS+GrLacb7jcfW0pbxO/LnunQJmkKIzLF3hZ4L4M4/8Ps7Ji26iqMdA5u4sSw4jBt3rXnH9x1WnVyVL/NuStAUQmRepebQYhyE/gohS0xa9Og27pQqYsWn64/zZuM3cSrqxHvb38t322NI0BRCZI3fWHBrDhvfNmn/ZnEbS8a2q07wpbvsPPWAj/w+IvBiIJv/Me3gU05J0BRCZI1ODz1mG1cOrRgIJtxhsmd9VzycizNl0yn6eQ6mcsnKvLftPQya6dPVZZcETSFE1hUvB91/gpvHjck9TESvU3zcuTbX7sUyb08Yk1pOIvRGKH1W9sk3GZAkaAohsse9jXEO56H5cGy1yYpt4FaKLl7lmbX7PMUty6FQLD+xPN+kjpOgKYTIvlYfgbMPbHgTHlw3WbFvBVQjPsnAzP0bUCgA4pLyR+o4CZpCiOzTW0L3WZAYC+vfMNn6dDeHorSs7sS5sIpYJaeO06HD383fJOXnhARNIUTOOLgbnzjPbILQpSYrdkATN2Kjq/Bx46V4OnlS1Koo9cvXN1n52SVBUwiRc41HgGtj2DTOuCWwCTSv6kBlx6IEnSrNZ60+417cPbad32aSsnNCgqYQIud0eug2HZLiYd1ok7ym63SKgU3cCA2LpKxVI0ralGTJMdNOqM9Wu/K6AUKIQqJ0FQj4FM5thSM/m6TIHvVcKGZtweIDV+lZqydrTq4h2gwJkbNCgqYQwnQaDDWuFto8Hh7cyHFxdtYW9Kzvysaj12hf5QWiE6JZf2a9CRqafRI0hRCmo9NBp2+Mo+lbPjRJkf19K5KkafwT7opzMWeW/J23r+gSNIUQpuVQFZq+AX8vhwu7c1ycm0NRWlV3YulfYfSq3ZtN5zZx++FtEzQ0eyRoCiFMr/lbUNLNmNQjMT7HxQ1o4satqHjKWwWQaEhk5YmVOW9jNknQFEKYnqUtdJgGt85A0Pc5Lq65uwNVHIsS+HdRajrUzNNRdAmaQgjzcA+Amp1h15dw91KOilLKOP3o2NX7tHDtzu5Luwm7F2aihmaNBE0hhPm0nwJKZ5JMSI+mHz282wiApcdMt/ooKyRoCiHMp4QL+L0NpzfCxb05KqqotQW9Griy97QF9co2YPbh2UzeMznXMx9J0BRCmFfjkVDc2TgFyZCzZMKPph/ZUYuzd87y4c4Pcz1lnARNIYR5WdpCqw/h6mE4nrO8mxVLF6W5uyPnbyUAYNAMxCfF52rKOAmaQgjzq9MbynjC9k8hMS5HRbWq7khidN2UPJtWeqtcTRknQVMIYX46PbSdCJGX4a/ZOSqqeTVHrA018XBojp2VHdte3oavq6+JGvp0EjSFELmjSiuo0hp2fwkxd7NdTGWHojjb21ICX6LioyhXrJwJG/l0EjSFELknYCLERkLQj9kuQilFc3cHrt2sAEBQuIyeCyEKq7IeUKsrHJgFD+9kuxi/ao4kxLlQxKIo+8L2mbCBTydBUwiRu1qMg7j7sH9GtotoUqU0eqXHuainPGkKIQq5MrWNyysPzMx236Z9ESvquNijEqoTej00VxMTS9AUQuS+lKfNmdkuws/dgXuRbiRpSRy8etCEjcuYBE0hRO4r6wk1Ohlf0WMis1WEXzVHLJNqAORqv6YETSFE3mjxLsTdg4PZm7fp5WqPvXVJSltXytV+TQmaQoi8Uc7LOG/zr9nZWiVkqdfhW6U0+oTqBIUFoZlgB8zMkKAphMg7vq9B1A04tipblzev5khirDu3Y25z9s5ZEzcubRI0hRB5p0orcKxpnOyejSfFFu6OWBtyt19TgqYQIu8oZXzavHEMLuzK8uUVShehasnqWGDL9IPTcyVFnARNIUTe8uwJRR2zvbSysvM1ErVYDl49mCu5NSVoCiHylqUNNBgKZ7dAxJksX65ZHQeMr/a5kVtTgqYQIu/VHww6Szg0P8uXdqkZAFgAYKm3NHtuTQmaQoi8Z+doXFoZsgQSYrN0afdaLSmrDQfg81afmz23pgRNIUT+4DPQmDbu5LosXabTKXycAgDjk6a55ShoKqUmKaWOKqVClFJblFLlk48rpdR3SqlzyZ/XS3XNAKXU2eSvAamO+yil/k6+5jullMpJ24QQBYxbcyhVGYKz/oru41IZpVlz9s45MzTscTl90vxS07Q6mqZ5AxuACcnHnwPck7+GATMAlFKlgI+BRkBD4GOlVMnka2Ykn/vouvY5bJsQoiDR6aDeALi8DyJOZ+nS2s4lsNDK8ff1rF2XHTkKmpqm3U/1bVEeDWFBV2CRZrQfsFdKlQPaAVs1TbujadpdYCvQPvmz4pqmBWnGtVCLgG45aZsQogDy7ps8ILQwS5fVLl8cS608Z2/n/ydNlFL/U0qFAX3590nTGQhLdVp48rGMjoencTy9OocppYKVUsERERE5vQUhRH5h5wg1OkDor5CUkOnL3J2KYUU5rkVfIsmQZMYGZiJoKqW2KaWOpfHVFUDTtPGaprkCi4HXH12WRlFaNo6nSdO0nzRNq69pWn1HR8en3YIQoiCp8yLE3IF/dmT6EisLHc52lUjSEgi/H/70C3LgqUFT07Q2mqZ5pPG19j+nLgGeT/5zOOCa6jMX4OpTjrukcVwI8ayp2gZsS8LR5Vm6rLZTNQDO3jZv4o6cjp67p/q2C3Aq+c/rgP7Jo+iNgXuapl0DNgNtlVIlkweA2gKbkz97oJRqnDxq3h/4b1AWQjwLLKygdnc4tRHiHmT6soauNQE4cs28g0E57dOckvyqfhRjAByTfPx34DxwDpgNjATQNO0OMAk4mPw1MfkYwAhgTvI1/wCbctg2IURB5dkLEmOMgTOTmlWuAZolweEnzNiwR2uPsknTtOfTOa4Br6Xz2TxgXhrHgwGPnLRHCFFIuDaCEhWMr+heL2bqktrl7bHUynLqVj5+PRdCCLPQ6aBOTzi/E6IyN0PGztqC4pauhD+4YN6mmbV0IYTIrtrdQTPA6cy/ope0teNu/AX2XTZfQmIJmkKI/KmMB9hXhJMbMnV6UFgQF6J3oJFE65/Nl1dTgqYQIn9Sypj56MIuiL331NMDLwZi0AwAxCeaL6+mBE0hRP5VszMkxcPZrU891d/NH0udMcuRTmdhtryaEjSFEPmXS0Mo6gQn1z/1VF9XX+Z0/hmAtq6DzZZXU4KmECL/0umgRkfjk2ZCzFNPf75WRwASEq3N1ySzlSyEEKZQoyMkRMOlvU89tYhVEfQU4Wb0TbM1R4KmECJ/c2sGFjZwdlumTrfVl+RO7C2zNUeCphAif7O0NQbOc5kLmsWsSnM//rbZmiNBUwiR/1VtA7fPwt2LTz21pI0jMYl3MBjSzS6ZIxI0hRD5X1XjxmmZedp0LOJIoorkVnScWZoiQVMIkf+VrmJcHZSJfs3yxcpg4D5XIx+apSkSNIUQ+Z9S4B4AF3Y/dRuMiiXLgzJwJuKaWZoiQVMIUTBU8jNOPbp6JMPTKpcqD8C5W+bZ/EGCphCiYKjYzPjfC7szPO1R0Lx4V540hRDPsqKlwak2XNyT4Wll7coAEHZPgqYQ4lnn1gwuH4DE+HRPcSrqBEDo7Q1mSQ8nQVMIUXBUam7cO+jKoXRPOX3buLHatfg/ab3I9Hk1JWgKIQqOik0BleEr+p5L/34Wn2T6vJoSNIUQBUeRUuBUE8IOpHuKMY+mAg2s9FYmz6spQVMIUbC4NIDwg2AwpPmxr6svZYu4YaG5sKTbRpPn1ZSgKYQoWFwbGre/uHUm3VMci5RBr9lTuURdk1cvQVMIUbC4NDT+N/yvdE8pYVMcTcVwJzr9UfbskqApdP7GCgAAEx1JREFUhChYSlcFG3sISz9o2tsWx0AMdx9K0BRCPOt0OuMrevjBdE9xsLXHoB5yJzrjderZqt7kJQohhLm5NICIUxB7P82PSxUpjkYMd+X1XAghgHLexv9eP5rmxyVsSqCpOG5F/3979x4j11necfz7m5m9ete3GMe3DbZTW60dijFLmi0FXKCJE7V1kPgjVCqmQk0VJZT+gUograAF1JKWIkUKkYziNlRVTUQpsVCiNARchGSSmMRJ7ATXGzvEjh3f75f1Xp7+cc5uZtczuzPeY3Zn9veRRmfmPWfOeV+fzZP3Pe85z2SfHs5B08xqz4I0aB7YXnJ1e2M7AIfPnMr80A6aZlZ72uZC+wI4+GLJ1e1NadA8dyLzQztomlltmv9uOFi6pzm9aToAx8+7p2lmlliwCo7uhp4zl60aHJ4fv1B6omg8HDTNrDbNXwUEvLXjslWDPc1TF08Tke2vUjpomlltmndDsjx0edAcvKbZF+c4fbEv08M6aJpZbZq+EJpmwOFXLl+V9jTP53/K06+Nnum9Wg6aZlabJLh2BRx+9bJVvzzySwDO57fyJz/4w0wTETtomlntmrsCDr0CI65bPnsgfS5dQe9AtomIHTTNrHZduwJ6TsHpN4cVf2TJR5I3IRpyDZkmInbQNLPaNXdlsjw0/Lrm+697P3kVaBpYyT984NFMExE7aJpZ7Zr7m8ny6K7LVrU0tNIYS1k6Y3Wmh3TQNLPa1TILpr2jZBb35nwzwSV6ekv/LMaVctA0s9o2Z3nyZNAILQ3NBL309PVnejgHTTOrbXOWle5pFtKeZp97mmZmb5uzHM4fg3PHhhW3FFoI9TpompkNc82yZHls+BA9GZ5foqfXw3Mzs7fN+Y1kOeK6ZnOhGeXc0zQzG27GdZArwIm9w4qbC83Iw3MzsxHyBZjRAcf3DCtuaWhJJ4I8PDczG272Ejh+eU8T+T5NM7PLzV5acnie3KfpoGlmNtysJXDxFJw/PlR08sJJeuMUb5x9IdNDOWiaWe2bvSRZpkP0rfu28sPdP6Sfizxx6G7n0zQzG2bmO5PlqTcA2PL6FvoHkgmg/uh1Pk0zs2FmdiTLk/sAWLN4DflcHoAczqdpZjZc84zk94JOJj3Nro4uPvXuTwHQ2foN59M0M7vMzA44tW/o49JZSwFoZlmmh3HQNLP6MPO6oeE5QCFXAOBi76VMD5NJ0JT0OUkhaU76WZIekNQt6SVJq4u2XS9pd/paX1T+Xkkvp995QJKyqJuZTREzhvc0G/INAFzsm2RBU1IH8AfAG0XFtwLL0tedwEPptrOBLwG/A9wIfEnSrPQ7D6XbDn5v7XjrZmZTyIxF0HM6uV8TaMglQbOnvzfTw2TR0/wm8NdA8W9orgO+E4mfAzMlzQduAZ6KiOMRcQJ4ClibrpseEVsjIoDvALdnUDczmyqmL0iWpw8Cb/c0e/p6Mj3MuIKmpD8G3oyIF0esWgjsK/q8Py0brXx/ifJyx71T0jZJ244cOTKOFphZ3RgKmsnP+Q72NHsH+ukfiHLfqlphrA0k/QiYV2LVfcAXgZtLfa1EWVxBeUkRsQHYANDZ2Zndv4aZ1a7BoHlmeE8z6ONS3wAtjflMDjNm0IyIj5Yql/QuYAnwYjpnswh4XtKNJD3FjqLNFwEH0vI1I8q3pOWLSmxvZlaZ9vnJ8nQSOgZ7mtBPT19/ZkHziofnEfFyRMyNiMURsZgk8K2OiLeAzcAn01n0m4BTEXEQeBK4WdKsdALoZuDJdN0ZSTels+afBB4bZ9vMbCopNEHrnLeDZlFPM8tMR2P2NK/Q48BtQDdwHvgzgIg4LukrwHPpdn8fEYNpSe4C/g1oAZ5IX2ZmlZs+fyhoDt6nGfRnmlMzs6CZ9jYH3wdwd5ntNgIbS5RvA27Iqj5mNgW1zYOzh4Di4Xlfptnb/USQmdWPtmvh7GGgaHiubIfnDppmVj/a5sK5wzAwMNTTPJv/Mc8d+Hlmh3DQNLP60XYtDPTBhePsOrYLgHP5p7jrydszS0TsoGlm9aNtbrI8e4jtb21P3ivoHcguEbGDppnVj/b0OZwzb9G5oDN5H6Ihl10iYgdNM6sf096RLM8d5T3z3pMU9X+I+z/0vcwSETtomln9mDYnWZ4/OvRzFy39XSyf9d7MDuGgaWb1o3km5Apw7ih5pY9NaiDThB0OmmZWP6TkUcpzR4Z6mjBA34Dv0zQzK23aHDh/bKinGbinaWZW3rRSPU0HTTOz0lrnDL+mSbZJiB00zay+tMyCCyeGeprhnqaZ2ShaZ8PFU+SH4uQA/f2eCDIzK61lNhDke86mBe5pmpmV1zobgPyl04Bnz83MRteSBs0Lp9MC9zTNzMprnQVAvucUABfyL9B98vnMdu+gaWb1pXkmANvS1HAXc7/gwR3rnU/TzKykNGj+bCifJvQ5n6aZWRnNMwD4UNuC5HOIgvNpmpmVkS9AYzu/25T0OKfFKj617GHn0zQzK6t5BrkLyUTQNN3AgtZ3Z7ZrB00zqz8tM9HFkwDkFL5P08xsVM0zUM8ZAKTwfZpmZqNqaoeeU+SUS3uafvbczKy8punQc4accki4p2lmNqrm6XDxdBo0fU3TzGx0TdOh5+2g6Z6mmdlomtphoG9oeN7f76BpZlZe83QABEjOcmRmNrrGdgByiJzw7LmZ2aia2gDISeBrmmZmY2hMgubAwAAne1/l0MWXMtu1g6aZ1Z+mNrbSx5m+8xzv28GPj37G+TTNzMpqbGML/UMf+3E+TTOz8hrbWEN+6GOe7PJpFjLZi5nZZNI4jS4KzCq0Ii3mXW2fzSyfpoOmmdWfxmnJQjnaG5ZzTcMNme3aw3Mzqz/5Bsg3IgCyu90IHDTNrF41TkME4aBpZlaBhmkowD1NM7NKpD3NrDlomll9amz18NzMrGINrSgCD8/NzCrR0JIOzx00zczG1tCCIoiML2s6aJpZfWpodU/TzKxiDS3kAk8EmZlVpKEVxQCQXdZ2cNA0s3pVaOYiA5zs3cOx3h2Z7dZB08zq0tZzhzgQwcm+bn524i+dhNjMbDRbzrwxdDXTSYjNzMawZs7KNMuRkxCbmY2pa+5KOhDn8tdxw/QvOgmxmdmoCs20IloKCyZPEmJJX5b0pqTt6eu2onVfkNQtaZekW4rK16Zl3ZLuLSpfIukZSbslfVdS43jqZmZTXKEZAE3CW46+GRGr0tfjAJJWAHcAK4G1wLck5SXlgQeBW4EVwCfSbQG+nu5rGXAC+HQGdTOzqarQjABF/5ibVuNqTQStAzZFRE9E7AW6gRvTV3dE7ImIS8AmYJ0kAR8Gvpd+/xHg9qtUNzObCgrJz13kJmFP8x5JL0naKGlWWrYQ2Fe0zf60rFz5NcDJiOgbUW5mdmUmangu6UeSdpR4rQMeAq4HVgEHgW8Mfq3EruIKysvV6U5J2yRtO3LkyFhNMLOpqNCUDs+zDZpjzp5HxEcr2ZGkbwM/TD/uBzqKVi8CDqTvS5UfBWZKKqS9zeLtS9VpA7ABoLOzM/t89mZW+/Jp0JxMw3NJ84s+fgwYfMBzM3CHpCZJS4BlwLPAc8CydKa8kWSyaHNEBPAT4OPp99cDj42nbmY2xRWagOyD5njv07xf0iqSofTrwF8ARMROSY8CrwB9wN0RyRSWpHuAJ4E8sDEidqb7+jywSdJXgReAh8dZNzObytLfPZ9UQTMi/nSUdV8Dvlai/HHg8RLle0hm183Mxi+9ppnL+Jqmnz03s/qUH3w+xkHTzGxshSbOE7w1cMz5NM3MxrL14PO8RnBg4ITzaZqZjWXLr35KpHeAO5+mmdkY1ixe43yaZmaV6uroYhl5zmsa18/8J+fTNDMbS5vyzFbb5MmnaWY2qSmH/LvnZmYVkibXs+dmZpOZkHuaZmYVk/ATQWZmFXJP08ysGhKEg6aZWUXk2XMzs2p4eG5mVjmJUX5u7Io4aJpZ3Uomgjx7bmZWGU8EmZlV7uxAH2/Q4yTEZmZj2bpvK7+8dIq99DoJsZnZWLa8vmXoaqaTEJuZjWHN4jVDAc5JiM3MxtDV0cVvNV/DhYun6Jj1gJMQm5mNpT3fxDzyzHASYjOzCshPBJmZVSV8n6aZWWWcsMPMrBqDw/MMe5sOmmZWx5JfPs9l+Py5g6aZ1b0CfZnty0HTzOqXkiuauXBP08xsTFIS4jw8NzOrQo7+DPdlZla3PBFkZlaxM/0X2csARy/tzGyfDppmVpe27tvKzjMHeI3gidN/43yaZmajSfJpBggG6HM+TTOz0ST5NAUBeQrOp2lmNpquji5Wtndw/vQ+VrX/rfNpmpmNZXpjK3PJsbDh+sz26eG5mdWx5JajvJ8IMjOrQBIzfXO7mVll0mfPfXO7mdnYpKSrmWUiYgdNM6tj6WOUvqZpZlY5X9M0M6uEnLDDzKwKnggyM6vY4ESQr2mamVVkcHju2XMzswoM3nKUXU/Tz56bWd063XuOMwywr29vZvt0T9PM6tLWfVt5+fhu9hD867lvOQmxmdlotry+JbmSKein30mIzcxGs2bxGprzTeQDmvLZJSF20DSzutTV0cXTf/QwX6GJpz/4d05CbGY2lq75q+miCWYvz2yf4+5pSvqMpF2Sdkq6v6j8C5K603W3FJWvTcu6Jd1bVL5E0jOSdkv6rqTG8dbNzKY4pSFustzcLun3gXXAb0fESuCf0/IVwB3ASmAt8C1JeUl54EHgVmAF8Il0W4CvA9+MiGXACeDT46mbmRm5fLKcLEETuAv4x4joAYiIw2n5OmBTRPRExF6gG7gxfXVHxJ6IuARsAtYpedbpw8D30u8/Atw+zrqZ2VSXPkY5mYLmcuAD6bD6fyW9Ly1fCOwr2m5/Wlau/BrgZET0jSg3M7tyV2F4PuZEkKQfAfNKrLov/f4s4CbgfcCjkpYy9MscwyTJRkqXl9u+XJ3uBO4EuO6660arvplNZYVmWNgJLbOz2+VYG0TER8utk3QX8P2ICOBZSQPAHJKeYkfRpouAA+n7UuVHgZmSCmlvs3j7UnXaAGwA6OzszO5JfDOrL+3z4M+fznSX4x2e/4DkWiSSlgONJAFwM3CHpCZJS4BlwLPAc8CydKa8kWSyaHMadH8CfDzd73rgsXHWzcwsc+O9T3MjsFHSDuASsD4NgDslPQq8AvQBd0dEP4Cke4AngTywMSJ2pvv6PLBJ0leBF4CHx1k3M7PMKYlxtauzszO2bds20dUwszoj6RcR0Tmy3I9RmplVwUHTzKwKDppmZlVw0DQzq4KDpplZFRw0zcyq4KBpZlYFB00zsyo4aJqZVcFB08ysCg6aZmZVcNA0M6uCg6aZWRUcNM3MquCgaWZWhZrPpynpCPCrCjadQ5JVvl7UU3vclslpqrflnRHxjpGFNR80KyVpW6mEorWqntrjtkxObktpHp6bmVXBQdPMrApTKWhumOgKZKye2uO2TE5uSwlT5pqmmVkWplJP08xs3KZE0JS0VtIuSd2S7p3o+lRL0uuSXpa0XdK2tGy2pKck7U6Xsya6nqVI2ijpsKQdRWUl667EA+l5eknS6omr+eXKtOXLkt5Mz812SbcVrftC2pZdkm6ZmFqXJqlD0k8kvSppp6TPpuU1d25GacvVOTcRUdcvIA+8BiwFGoEXgRUTXa8q2/A6MGdE2f3Aven7e4GvT3Q9y9T9g8BqYMdYdQduA54ABNwEPDPR9a+gLV8GPldi2xXp31oTsCT9G8xPdBuK6jcfWJ2+bwf+L61zzZ2bUdpyVc7NVOhp3gh0R8SeiLgEbALWTXCdsrAOeCR9/whw+wTWpayI+ClwfERxubqvA74TiZ8DMyXN//XUdGxl2lLOOmBTRPRExF6gm+RvcVKIiIMR8Xz6/gzwKrCQGjw3o7SlnHGdm6kQNBcC+4o+72f0f9DJKID/kfQLSXemZddGxEFI/miAuRNWu+qVq3utnqt70iHrxqLLJDXTFkmLgfcAz1Dj52ZEW+AqnJupEDRVoqzWbhl4f0SsBm4F7pb0wYmu0FVSi+fqIeB6YBVwEPhGWl4TbZHUBvwX8FcRcXq0TUuUTar2lGjLVTk3UyFo7gc6ij4vAg5MUF2uSEQcSJeHgf8mGUocGhwepcvDE1fDqpWre82dq4g4FBH9ETEAfJu3h3mTvi2SGkiCzH9ExPfT4po8N6XacrXOzVQIms8ByyQtkdQI3AFsnuA6VUzSNEntg++Bm4EdJG1Yn262HnhsYmp4RcrVfTPwyXSm9ibg1OBQcbIacV3vYyTnBpK23CGpSdISYBnw7K+7fuVIEvAw8GpE/EvRqpo7N+XactXOzUTPfP2aZtduI5lRew24b6LrU2Xdl5LM9L0I7BysP3AN8DSwO13Onui6lqn/f5IMjXpJ/g//6XJ1Jxk2PZiep5eBzomufwVt+fe0ri+l/zHOL9r+vrQtu4BbJ7r+I9ryeyRD0peA7enrtlo8N6O05aqcGz8RZGZWhakwPDczy4yDpplZFRw0zcyq4KBpZlYFB00zsyo4aJqZVcFB08ysCg6aZmZV+H+ZoBDzLsOchgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,8))\n",
    "ax.plot(dz_new, -eh_new, label='original')\n",
    "ax.plot(ds.dz, -ds.z, label = 'initial guess using func_tanh')\n",
    "ax.plot(func_tanh(eh_new, *popt), -eh_new, 'g.-', label='best fit using curve_fit')\n",
    "ax.legend()\n",
    "#ax.set_ylim(-100,0)\n",
    "#ax.set_xlim(0,30)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.94,   2.09,   2.24,   2.39,   2.77,   3.38,   4.01,   4.65,\n",
       "         5.29,   5.95,   6.61,   7.28,   7.97,   8.66,   9.37,  10.08,\n",
       "        10.81,  11.54,  12.29,  13.06,  13.85,  14.69,  15.59,  16.56,\n",
       "        17.61,  18.76,  20.02,  21.42,  23.  ,  24.77,  26.79,  29.1 ,\n",
       "        31.76,  34.87,  38.5 ,  42.79,  47.9 ,  54.01,  61.37,  70.25,\n",
       "        80.95,  93.75, 108.8 , 126.04, 145.04, 164.81, 184.05, 201.34,\n",
       "       215.66, 226.64, 234.5 , 239.84, 243.31, 245.52, 246.88, 247.72,\n",
       "       248.23, 248.54, 248.73, 248.84, 248.91, 248.95, 248.98, 248.99,\n",
       "       249.  ])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's truncate the array using 2 decimal digits\n",
    "dz_best_fit = func_tanh(eh_new, *popt).round(decimals=2)\n",
    "dz_best_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6000.01"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_best_fit.cumsum()[-1].round(decimals=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) let assure the sum of the first four levels equals MINIMUM DEPTH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  2.5 ,   2.5 ,   2.5 ,   2.5 ,   2.77,   3.38,   4.01,   4.65,\n",
       "         5.29,   5.95,   6.61,   7.28,   7.97,   8.66,   9.37,  10.08,\n",
       "        10.81,  11.54,  12.29,  13.06,  13.85,  14.69,  15.59,  16.56,\n",
       "        17.61,  18.76,  20.02,  21.42,  23.  ,  24.77,  26.79,  29.1 ,\n",
       "        31.76,  34.87,  38.5 ,  42.79,  47.9 ,  54.01,  61.37,  70.25,\n",
       "        80.95,  93.75, 108.8 , 126.04, 145.04, 164.81, 184.05, 201.34,\n",
       "       215.66, 226.64, 234.5 , 239.84, 243.31, 245.52, 246.88, 247.72,\n",
       "       248.23, 248.54, 248.73, 248.84, 248.91, 248.95, 248.98, 248.99,\n",
       "       249.  ])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_best_fit[0:4] = 2.5\n",
    "dz_best_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6001.35"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_best_fit.cumsum()[-1].round(decimals=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) let assure that $\\sum{dz}$ = MAXIMUM DEPT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.35"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_offset = (dz_best_fit.cumsum()[-1] - h_max).round(decimals=2)\n",
    "dz_offset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's deflate the bottom-most five levels by dz_offset/5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "dz_best_fit[-5::] = (dz_best_fit[-5::] - dz_offset/5.).round(decimals=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6000.0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_best_fit.cumsum()[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  2.5 ,   2.5 ,   2.5 ,   2.5 ,   2.77,   3.38,   4.01,   4.65,\n",
       "         5.29,   5.95,   6.61,   7.28,   7.97,   8.66,   9.37,  10.08,\n",
       "        10.81,  11.54,  12.29,  13.06,  13.85,  14.69,  15.59,  16.56,\n",
       "        17.61,  18.76,  20.02,  21.42,  23.  ,  24.77,  26.79,  29.1 ,\n",
       "        31.76,  34.87,  38.5 ,  42.79,  47.9 ,  54.01,  61.37,  70.25,\n",
       "        80.95,  93.75, 108.8 , 126.04, 145.04, 164.81, 184.05, 201.34,\n",
       "       215.66, 226.64, 234.5 , 239.84, 243.31, 245.52, 246.88, 247.72,\n",
       "       248.23, 248.54, 248.73, 248.84, 248.64, 248.68, 248.71, 248.72,\n",
       "       248.73])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz_best_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare old and new grids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# new interfaces\n",
    "e = np.zeros(len(dz_best_fit)+1)\n",
    "e[1::] = dz_best_fit.cumsum()\n",
    "eh_best_fit = (0.5 * (e[0:-1]+e[1::])).round(decimals=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><table><tr><th>Level #</th><th>new layer depth (m)</th><th>new thickness (m)</th><th>old layer depth (m)</th><th>old thickness (m)</th></tr><tr><td>1</td><td>1.25</td><td>2.5</td><td>1.25</td><td>2.5</td></tr><tr><td>2</td><td>3.75</td><td>2.5</td><td>3.75</td><td>2.5</td></tr><tr><td>3</td><td>6.25</td><td>2.5</td><td>6.25</td><td>2.5</td></tr><tr><td>4</td><td>8.75</td><td>2.5</td><td>8.75</td><td>2.5</td></tr><tr><td>5</td><td>11.385</td><td>2.77</td><td>15.0</td><td>10.0</td></tr><tr><td>6</td><td>14.46</td><td>3.38</td><td>25.0</td><td>10.0</td></tr><tr><td>7</td><td>18.155</td><td>4.01</td><td>35.0</td><td>10.0</td></tr><tr><td>8</td><td>22.485</td><td>4.65</td><td>45.0</td><td>10.0</td></tr><tr><td>9</td><td>27.455</td><td>5.29</td><td>55.0</td><td>10.0</td></tr><tr><td>10</td><td>33.075</td><td>5.95</td><td>65.0</td><td>10.0</td></tr><tr><td>11</td><td>39.355</td><td>6.61</td><td>75.0</td><td>10.0</td></tr><tr><td>12</td><td>46.3</td><td>7.28</td><td>85.0</td><td>10.0</td></tr><tr><td>13</td><td>53.925</td><td>7.97</td><td>95.0</td><td>10.0</td></tr><tr><td>14</td><td>62.24</td><td>8.66</td><td>105.0</td><td>10.0</td></tr><tr><td>15</td><td>71.255</td><td>9.37</td><td>115.0</td><td>10.0</td></tr><tr><td>16</td><td>80.98</td><td>10.08</td><td>125.0</td><td>10.0</td></tr><tr><td>17</td><td>91.425</td><td>10.81</td><td>135.0</td><td>10.0</td></tr><tr><td>18</td><td>102.6</td><td>11.54</td><td>145.0</td><td>10.0</td></tr><tr><td>19</td><td>114.515</td><td>12.29</td><td>155.0</td><td>10.0</td></tr><tr><td>20</td><td>127.19</td><td>13.06</td><td>165.125</td><td>10.25</td></tr><tr><td>21</td><td>140.645</td><td>13.85</td><td>175.5</td><td>10.5</td></tr><tr><td>22</td><td>154.915</td><td>14.69</td><td>186.25</td><td>11.0</td></tr><tr><td>23</td><td>170.055</td><td>15.59</td><td>197.625</td><td>11.75</td></tr><tr><td>24</td><td>186.13</td><td>16.56</td><td>209.75</td><td>12.5</td></tr><tr><td>25</td><td>203.215</td><td>17.61</td><td>222.625</td><td>13.25</td></tr><tr><td>26</td><td>221.4</td><td>18.76</td><td>236.375</td><td>14.25</td></tr><tr><td>27</td><td>240.79</td><td>20.02</td><td>251.25</td><td>15.5</td></tr><tr><td>28</td><td>261.51</td><td>21.42</td><td>267.5</td><td>17.0</td></tr><tr><td>29</td><td>283.72</td><td>23.0</td><td>285.375</td><td>18.75</td></tr><tr><td>30</td><td>307.605</td><td>24.77</td><td>305.0</td><td>20.5</td></tr><tr><td>31</td><td>333.385</td><td>26.79</td><td>326.75</td><td>23.0</td></tr><tr><td>32</td><td>361.33</td><td>29.1</td><td>351.0</td><td>25.5</td></tr><tr><td>33</td><td>391.76</td><td>31.76</td><td>378.125</td><td>28.75</td></tr><tr><td>34</td><td>425.075</td><td>34.87</td><td>408.75</td><td>32.5</td></tr><tr><td>35</td><td>461.76</td><td>38.5</td><td>443.375</td><td>36.75</td></tr><tr><td>36</td><td>502.405</td><td>42.79</td><td>482.75</td><td>42.0</td></tr><tr><td>37</td><td>547.75</td><td>47.9</td><td>527.75</td><td>48.0</td></tr><tr><td>38</td><td>598.705</td><td>54.01</td><td>579.375</td><td>55.25</td></tr><tr><td>39</td><td>656.395</td><td>61.37</td><td>638.875</td><td>63.75</td></tr><tr><td>40</td><td>722.205</td><td>70.25</td><td>707.625</td><td>73.75</td></tr><tr><td>41</td><td>797.805</td><td>80.95</td><td>787.125</td><td>85.25</td></tr><tr><td>42</td><td>885.155</td><td>93.75</td><td>879.0</td><td>98.5</td></tr><tr><td>43</td><td>986.43</td><td>108.8</td><td>984.875</td><td>113.25</td></tr><tr><td>44</td><td>1103.85</td><td>126.04</td><td>1106.375</td><td>129.75</td></tr><tr><td>45</td><td>1239.39</td><td>145.04</td><td>1244.75</td><td>147.0</td></tr><tr><td>46</td><td>1394.315</td><td>164.81</td><td>1400.625</td><td>164.75</td></tr><tr><td>47</td><td>1568.745</td><td>184.05</td><td>1574.0</td><td>182.0</td></tr><tr><td>48</td><td>1761.44</td><td>201.34</td><td>1764.0</td><td>198.0</td></tr><tr><td>49</td><td>1969.94</td><td>215.66</td><td>1968.875</td><td>211.75</td></tr><tr><td>50</td><td>2191.09</td><td>226.64</td><td>2186.375</td><td>223.25</td></tr><tr><td>51</td><td>2421.66</td><td>234.5</td><td>2413.875</td><td>231.75</td></tr><tr><td>52</td><td>2658.83</td><td>239.84</td><td>2648.875</td><td>238.25</td></tr><tr><td>53</td><td>2900.405</td><td>243.31</td><td>2889.25</td><td>242.5</td></tr><tr><td>54</td><td>3144.82</td><td>245.52</td><td>3133.25</td><td>245.5</td></tr><tr><td>55</td><td>3391.02</td><td>246.88</td><td>3379.625</td><td>247.25</td></tr><tr><td>56</td><td>3638.32</td><td>247.72</td><td>3627.5</td><td>248.5</td></tr><tr><td>57</td><td>3886.295</td><td>248.23</td><td>3876.25</td><td>249.0</td></tr><tr><td>58</td><td>4134.68</td><td>248.54</td><td>4125.5</td><td>249.5</td></tr><tr><td>59</td><td>4383.315</td><td>248.73</td><td>4375.125</td><td>249.75</td></tr><tr><td>60</td><td>4632.1</td><td>248.84</td><td>4625.0</td><td>250.0</td></tr><tr><td>61</td><td>4880.84</td><td>248.64</td><td>4875.0</td><td>250.0</td></tr><tr><td>62</td><td>5129.5</td><td>248.68</td><td>5125.0</td><td>250.0</td></tr><tr><td>63</td><td>5378.195</td><td>248.71</td><td>5375.0</td><td>250.0</td></tr><tr><td>64</td><td>5626.91</td><td>248.72</td><td>none</td><td>none</td></tr><tr><td>65</td><td>5875.635</td><td>248.73</td><td>none</td><td>none</td></tr><tr><td>Total thickness</td><td></td><td>6000.0</td><td></td><td>5500.0</td></table></html>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import HTML, display\n",
    "\n",
    "html = \"\"\"<html><table><tr><th>Level #</th><th>new layer depth (m)</th><th>new thickness (m)</th><th>old layer depth (m)</th><th>old thickness (m)</th></tr>\"\"\"\n",
    "for k in range(len(dz_new)):\n",
    "    html += \"<tr><td>{}</td>\".format(str(k+1))\n",
    "    html += \"<td>{}</td>\".format(str(eh_best_fit[k]))\n",
    "    html += \"<td>{}</td>\".format(str(dz_best_fit[k]))\n",
    "    if k < len(old_vgrid.dz):\n",
    "        html += \"<td>{}</td>\".format(str(eh[k]))\n",
    "        html += \"<td>{}</td>\".format(str(old_vgrid.dz[k].values))\n",
    "    else:\n",
    "        html += \"<td>none</td>\"\n",
    "        html += \"<td>none</td>\"\n",
    "        \n",
    "    html += \"</tr>\"\n",
    "html += \"<tr><td>Total thickness</td><td></td><td>{}</td><td></td><td>{}</td>\".format(dz_best_fit.cumsum()[-1],old_vgrid.dz.cumsum()[-1].values)    \n",
    "html += \"</table></html>\"\n",
    "display(HTML(html))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save new vertical grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.Dataset(data_vars={ 'dz' : (('z'), dz_best_fit)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "ds.dz.attrs['long_name'] = 'Nominal thickness of level'\n",
    "ds.dz.attrs['units'] = 'meter'\n",
    "# Global attrs\n",
    "ds.attrs['title'] = 'Vertical resolution using {} layers'.format(len(dz_best_fit))\n",
    "ds.attrs['maximum_depth'] = '{} m'.format(dz_best_fit.cumsum()[-1])\n",
    "ds.attrs['minimum_depth'] = '10 m'\n",
    "ds.attrs['url'] = 'created using https://gist.github.com/gustavo-marques/3e11fd764749599559c45807c3f9b570'\n",
    "ds.attrs['author'] = 'Gustavo Marques (gmarques@ucar.edu)'\n",
    "ds.attrs['Created'] = datetime.now().isoformat()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = 'vgrid-{}L-{}.nc'.format(len(dz_best_fit), datetime.now().isoformat()[0:10])\n",
    "ds.to_netcdf(fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "analysis",
   "language": "python",
   "name": "analysis"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}