PFIO_tools_compare.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import time\n",
    "import pf_pytools.pf_fort_io as pf_fort_io\n",
    "import pfio\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_pfb = '/home/arezaii/git/tools/tests/press.init.pfb'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "load time 0.0027 seconds\n"
     ]
    }
   ],
   "source": [
    "t2 = time.perf_counter()\n",
    "pfb_c_data = pfio.pfread(input_pfb)\n",
    "# Flip not needed with new PFIO code\n",
    "pfb_c_data = np.flip(pfb_c_data,axis=1)\n",
    "t3 = time.perf_counter()\n",
    "print(f'load time {t3-t2:0.4f} seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "load time 0.0009 seconds\n"
     ]
    }
   ],
   "source": [
    "t2 = time.perf_counter()\n",
    "pfb_f_data = np.zeros((41,41,50), order='F')\n",
    "pf_fort_io.pfb_read(pfb_f_data, input_pfb)\n",
    "# X Y Z -> Z Y X\n",
    "pfb_f_data = np.transpose(pfb_f_data,(2,1,0))\n",
    "t3 = time.perf_counter()\n",
    "print(f'load time {t3-t2:0.4f} seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use numpy testing to make sure the data are all equal.\n",
    "\n",
    "this function returns None if the arrays are exactly the same, otherwise an exception will be thrown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_array_equal(pfb_f_data, pfb_c_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare for plotting\n",
    "\n",
    "Used xarray because it provides convenient plotting features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "xds1 = xr.DataArray(pfb_c_data)\n",
    "xds2 = xr.DataArray(pfb_f_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f2eae2d65d0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAEmCAYAAABRZIXOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZQc5Xkv/u/TPT2bZkajmdEyaEECJEBswsgyhIAxeAFiszh2rnHi4IQb2b8TJ3ac3BgnObGdxDn2vXbwPTe/OJEvNtg3Pww29oGLgzEBDDExwkIItCGEFrSNds1Io9l6up/fH90SM/O8Ule3qqbf6vp+zumD+pmarqruni9vV9dTr6gqiIiIiChcqWpvABEREVEt4iCLiIiIKAIcZBERERFFgIMsIiIioghwkEVEREQUAQ6yiIiIiCIwKYMsEUmLyMsi8ljx/gIRWSkim0XkQRGpn4ztICIqF/OLiColk3GdLBH5LIClANpU9f0i8hCAH6nq90XknwG8oqrfPN1j1E9t0sZZbSfvL2rsq3h7XHu8L9dgau2pQVNrlHDHpaPIm9qebKupNadGTK0rna1onTnHOl873mVq0pc2tbx9mrCoY6+p1Yv9XdfzfjRvn88Rtb/bIKOm1pqy+yEQx3rDfo/bdfTm7TanHc9zayrItthlsmofawR2nQCwae3wQVWdHmBFAID3vWuKHjqcC7o4AOClV4efUNUby/qlmIoiv4BwMyxJ+QX4k2HMr1OpPMPKzS+g/AybzPyqi3oFIjIHwG8A+DKAz4qIALgewEeLi9wP4IsAThtSjbPa8I5/+ejJ+09c8JOKtynvePP8j8OLTO3W1ldM7YLMlIrX63Iof9zUvrD3Xab2tpY3Te332/ZVtM6jeRu+v/biXabW+NOp9nfPtX88j97xP0xtXp0NWtfz/vhAs6ntyU4ztXPq95vatY3DppYR+5bOqg24M5F2/I/q4X67ze1p+9re0GT/ZzNRzhFGPTn7mu0cdb8Xr12wxb5ZTuPg4RxWPjGnnF9BpnuL/T9aDYoqv4BwMyxJ+QX4k2HML7czybBy8wsoP8MmM78iH2QB+AaAPwdw4h3bCaBX9eQ7ZxeA2ZOwHUTkpM5QJADML6IY8DfDIj0nS0TeD2C/qr40tuxY1HkMUkSWi8gqEVmV7bOjYCI6cwogDy3rlgTML6J4KDfDJlPUR7KuBnCLiNwMoBFAGwqfDNtFpK74aXAOgD2uX1bVFQBWAEDb+TOTkexEVeD6KpeYX0Rx4WuGRXokS1U/r6pzVHU+gI8AeFpVfxvAMwA+VFzsTgCPRLkdRHRqCkVOy7slAfOLKB7KzbDJNBnnZLl8DsD3ReTvALwM4N5SvyAC1I3pxvinXnuS27OHzze1EUfXxMBoxtRu715jaosyTabmGi2/e/1tprbtzRmmlm6yJzBmGmzN9X3E4sX2w7JrW1IBxs3PDnWaWs7RJTNsz4XEaLN9g97Xu8zU5tUfMrW+nD1J9DtvXGlqN87baGofaNlsammxj+d6Tl7L2i6m3rx9ba9uDPZJyPXd/z/vuNbU3j1zk6nd0PRaycd3nZg6p851wrI9MbVSSfkKMCRnnF9AuBmWpPwC/Mkw5pcbM+wtkzbIUtWfA/h58d9bAdh3NRFNOgWQ8zSgfMH8IvKXzxlWrSNZROQRXz8FEhEF4WuGcZBFlHAKJOY8KyKqPT5nGAdZRORpXw4RUTC+ZhgHWUQJp1Bvz2cgIirF5wyLzSArIznMaOw/ef9/bbjOLDPcY7sXGg/YLoehc+x0BnPn2U4Sl28cOdfUdh1qN7XmbbaDMZW1NVcrTt3Vh03tmmbbnQLYSbhcUzCMYvycThdk+s0ybU1Dpnao29EN0mznh3q5d66pXX6WnRlhSO3+33z2BlN7V6utzUq3mJqrE+cnA3a51pTtzrq6wTVvmn2v7Bi13S//88B1prZtt51qa+dUO5WGXaNdZ9Drvbi7dSqgQM7PfKoZE/MLCDfDkpRfgD8ZlvT8KqyVGXY6sRlkEVE0CldLJiKKJ58zjIMsosQT5JxXNyIiigN/M4yDLKKEUwB5Tw+1ExGV4nOGcZBFRN5+CiQiCsLXDOMgiyjhCldL9jOgiIhK8TnDYjPIakqPYEnrzpP3V+Jss8yUXbbLIWUbCVE/baDi7VjcuNvUzp5uu2n2iO2aaOizj3f4HSOmdkWH7RSamrJdNynY+atc77M/73nHuPuZlD1F8Ohgo6k1nX3M1P72YjsX7i1T7HLujhPbEQTsdNQs10mNG7O2m2jT0EJTe1eLnUssLcHe+jtHW01t/7CtpTN2C4+M2PfAsI7vCmoS22EVVNB53oLIq58BVSsm5hdQnQyrhfwCfMqwZOcXwAwrJTaDLCKKhs+fAomISvE5w8IbRhJRLCkEOaTKupUiIp8WkXUisl5EPlOsfVFEdovImuLt5sh3johqXrkZFkRYGcYjWUQU6qF2EbkYwB8AWAZgBMBPReQnxR/fo6pfC21lRETwN8M4yCJKuAgOtV8I4AVVHQAAEXkWwO1hroCI6ASfM4xfFxIlnEKQ1bqybiWsA3CtiHSKSDOAmwGcmLfkUyLyqoh8W0SmRblfRJQM5WZYAKFlWGyOZB3PNeCF3nNO3m9rtp0ZhzvsvE+prB3dXjXTdtjc1Gy7dVydD+91zI/14vQtpvbQqJ0LK+uYpmnGLNuyc0On7SZxzfEUdH6o+Y3ju32+ueEas8zIYduZM+e8HlN7X7PtREo55iALum27HHNrNYp9zTZnm01tz2i3qW0f6jS19jbbAZVCvam5tvnyevt6XzPNzsM2NTNoauc37zO1rI5fR9MZfPgK+hwHUcGnwC4RWTXm/gpVXQEAqrpRRL4K4EkA/QBeATAK4JsA/haFD55/C+DrAH7/DDc9FibmFxBuhiUpvwB/Mizp+QXENsNOmV9AuBkWm0EWEUVDVZDTsg9qH1TVpad+TL0XwL0AICJ/D2CXqp5MbRH5FoDHKthcIqJxKsiw0+ZX4THDyTAOsogI+ZDbn0VkhqruF5F5AD4I4CoR6VbVE4cWbkfhkDwR0RnzNcM4yCJKuMJJo6GfnvmwiHQCyAL4Q1U9IiLfE5ElxVVuB/CJsFdKRMnjc4ZxkEWUeBV9XXhaqmpOmlHVj4W6EiIiAD5nWKSDLBFpBPAcgIbiun6oql8QkfsAvBPAibMmP66qa6LcFiJyUwB5NhobzC+iePA5w6I+kjUM4HpV7ReRDIBfiMjjxZ/9N1X9YdAHGjjWiJd+sejk/VyTmmWaeu13skMX266J/zrj2UDr3D5q57Pakm03NdcIerjLbl9q2G7fpdMOmNrNU15zbI3tznF1tsyrs/NSXdw4fn6tfM5uh6bs9jbX2XmqNmXtcs8en29qvzhynqld0LLX1B5+Y4mpNdTb9V4/x9ERU2df2xn19jVbUGc7e1y2ZG2H1sx02tQW1tv9uKjTdqwurrevz8R5vlzdNTlHB8+hvO0Syji6mCqV83TeryqLLL+A6DOsVvML8CfDkp5fADOslEgHWaqqwMkZNTPFm32HE1HVnJiSgsZjfhHFg88ZFvlWiUhaRNYA2A/gSVVdWfzRl4sX9LpH5Aym8SaiM5bXVFm3pGB+EcWDr/kV+dpUNaeqSwDMAbCsOCfQ5wFcAODtADoAfM71uyKyXERWiciq3HF76JKIztyJzpywJ1itBcwvIv+Vm2GTadLWpqq9AH4O4EZV7dGCYQDfQWESRtfvrFDVpaq6ND3FcblhIjpjCkFOy7slDfOLyF/lZthkinSQJSLTRaS9+O8mAO8G8JqIdBdrAuA28KKERFWVR6qsWxIwv4jiw9f8irq7sBvA/SKSRmFA95CqPiYiT4vIdAACYA2AT5Z6oLrjwKwX3upYyDXY0Wiu3p6T2nat7da4pH7Y1F4atk/FP+272dR6R5pMbf+AnW8sPWS3b+QsOwfV78543tRc83y5TE/b+atc3tU0vtvl60ttU9T399sP45sOzzC1v9xuJyIfzGZMbc+Rqab2WsY+3uARO9/Y8IDd/0eOXWpqdfU5U/vgwldM7YFjXaY2o+6oqeXQYWrPD9purLkZO5fa1Y2uObjse6U/P/69t2nEvu/WD882tX2jbaZ2KGs7sQoeOkXdTdXdYUbR5RcQboYlKb8AfzIs6fkFRJFh5eUX4HeGRd1d+CqAyx3166NcLxGVQ0KfkqIWML+I4sLfDOMV34kSTuHvp0AiolJ8zjAOsogoUR2DRFR7fM0wDrKIEk4hyCewY5CIaoPPGcZBFhF5+ymQiCgIXzMsNoMsrQOGpr31JIptzMBImx3JtqRcXRPWS0PzTW1Vz1xTGzhmL+4sB22tzjH5RnO7natqVtp2DhXmoy3tcN52+2wYsZ0t52R6x91/eeAis8wre88ytaEB2/0zrcnOj3XHnF+ZWnqufd5XHZtvauct2m9qrx+fZWo/W7fY1EZNBbioaZepPX90oanNarDdOX/Vtd7UHle73GxHZ09f3r73pqZs51FLavxr25qy74nOun5Ty4jd28uadpgaAHzdWT01BRJ1FfdqmJhfQLgZlqT8AvzJsKTnFxB+hpWbX4DfGRabQRYRRUWQ87Qzh4ioNH8zjIMsooTz+VMgEVEpPmcYB1lE5O2nQCKiIHzNMA6yiBJOVbz9FEhEVIrPGcZBFhF5eyE/IqIgfM0wDrKIEk4Bb6ekICIqxecMi80gS9tzGL39yMn7A4O2NbexwU4kWp+2baO3v/ZbjuVsP/XxPjtBZscv7XpTWdvv3G/nx0RHi20dbk25Gnltm2zKcQ2QvTm73IahOaY2PT2+nXbjsW6zTH6NnQw1P8M+J3WzbVtzX67Z1P6sY4upfbytx9Scptnf/caUvab2+F7bxn1Bg13ujtkHTe2+o/Y52DV63NRusruGI46O+g0jtt15WYNdMCPj/+TOrbOv4aKMbX/OOVqx0xLWJzfx9lNgrZiYX0C4GZak/AI8z7AE5RfADCslNoMsIoqGAshqutqbQURUEZ8zjIMsooTzeUoKIqJSfM4wDrKICHlPp6QgIgrC1wzjIIso4VSBnKefAomISvE5wzjIIiJvD7UTEQXha4bFZpC1oOkgvnfpd07ezzo6CVydLnfvuM3U1j9pJ9x0zF2JujbbdZO1DTuok2Av7gXtdjLRmWnb7eMyqMOmtjPbYWov9J5jartHpo2731xnH2togZ2sNXUkY2qv75xpah0Ntuvoe5k+U/t57/mmdl37JlP7WKvtpvmMo2Pno21rTW39iO0wAmzHlqtL6PmhVlObU2c7bKal7Jvgsnr7/KWl9GsbtLvG1dUTlsL5DH4eaq8VE/MLCDfDkpRfgD8ZlvT8KizHDDud2AyyiCg6vk5JQUQUhK8ZxkEWUcIVJlf1M6CIiErxOcM4yCJKPH8PtRMRleZvhkW6VSLSKCIvisgrIrJeRL5UrC8QkZUisllEHhQJ+OUvEUUiDynrVoqIfFpE1hX/7j9TrHWIyJPFv/snRcSebOMR5hdRfISZX0B4GRb10G8YwPWqehmAJQBuFJErAXwVwD2quhDAEQB3RbwdRHQKJ9qfy7mdjohcDOAPACwDcBmA94vIQgB3A3iq+Hf/VPG+z5hfRDFQboaVEmaGRfp1oaoqgBOTGGWKNwVwPYCPFuv3A/gigG+e7rEaJYULM3Z+pbF6crYLY2jU7mK2xXbdpLL2iZf5tuMkvdh2tRw50GI3JmvHr0//x6WmdlHHhfZ37eYh3Wjn4EpnbC2/xW7LylnjO0c6uo7ZFQzaKQmm7LLPyXGx81Ttn23X+e+HF5va7KZeU3N1HKZgO2LucHTszErb9XY22tfnmUG7zTc02efu6kbHpF4BNaccc8KF+Bkm73hOwhTyofYLAbygqgMAICLPArgdwK0Arisucz+AnwP4XJgrDtNk5xdQeYYlKb8AfzKM+RVcUjMs8i8xRSQtImsA7AfwJIAtAHpV9UTT8S4AjulIiWgynJiSopwbgC4RWTXmtnzMQ64DcK2IdIpIM4CbAcwFMFNVewCg+N8Zk72v5WJ+Efmv3AzD6fMLCDHDIj/xXVVzAJaISDuAH6MwQjSLuX63uOPLAWDebJ6jTxSVoOcpjHFQVZe6fqCqG0XkqygMSvoBvALAcSUn/zG/iOKhzAw7ZX4B4WbYpJ2Or6q9KBxauxJAu8jJK5PNAbDnFL+zQlWXqurSrk4/Z9gmirsT7c9lfhI8/WOq3quqb1PVawEcBrAZwD4R6QaA4n/t1S09xfwi8le5GRboMUPKsKi7C6cXPwFCRJoAvBvARgDPAPhQcbE7ATwS5XYQ0enlNVXWrRQRmVH87zwAHwTwAIBHUfh7B2Lwd8/8IoqPMPMLCC/Doj6G3Q3gfhFJozCge0hVHxORDQC+LyJ/B+BlAPdGvB1EdCplfLorw8Mi0onCnCB/qKpHROQrAB4SkbsA7ADw4bBXGjLmF1EceJxhUXcXvgrgckd9KwqtkaH6z8E5ptY7bOdpynXYr1YXnbvb1D4792emtmHYnuO6cVa3qT3zpp0fMb3dzi1Vv9V2dTTvs10YmUE7+h6aZruVUqOuDo7xL/Nog50zbJZtagLUPtax823tmul2Xq4ZmaOmdmWTXW562q54xeErTe1/jbaZ2qenbTc11/xYTxy1XVGz6p43tYsyzaYWtokdNn35IbPMpqx9T6wdmmtq7WnbPVYJRUXnZJ3+MVWvcdQOAbgh1BVFaLLzC6g8w5KUX4A/GZb0/AKYYaXwbEwi8nZKCiKiIHzNMA6yiBLO53m/iIhK8TnDOMgiIm8DiogoCF8zjIMsooQ7cSE/IqI48jnDOMgiotBPGiUimky+ZlhNDbLWOTpzeg5MNbX6Vjs/1G/OWm1qF9TbufZeGZpnanMbj5jaeV12rqotV5gSBm0Jx1+3XTzNPcGu7ZGdYt9oZz0/fi2ZNw+YZQYW2w6jkan2AooNe+1bJue47sh1za+b2qKM7ZICbBfKF6avNbXXs/aZyqqd0ystdluub1tvakHmkStHpfNyTU3Z7WhP2X2dW3/I1I7lXM9nBdTfQ+1JU2mGJSm/AJ8yLNn5BTDDSqmpQRYRlc/nk0aJiErxOcM4yCIibwOKiCgIXzOMgyyihPP5pFEiolJ8zjAOsogI6mlAEREF4WuGcZBFlHCqwGjASVOJiHzjc4bFZpB1XBUvDr/VAdGbt/M09QzbLhxHswYymZypzazrM7VGxy+/c8omU9s6Mt3Ups6wczLtm2a3b93Rs0xtzWDG1Ppm2i4ZV5dktsc+L/m68fsxOrvTLNO7yK6z9xI7L1f7TNuxtOGY7erZ3WLn6lqUsY+XQrA/jAsyU0wtaEfMe5vs3FpB1+tax0B+xNRW9F1oaosb7XxyNzbZ353o3DrbdXQoZ/dhVr2dW61Svn4KrBUT8wsIN8OSlF+APxnG/HJjhr0lNoMsIoqKv+czEBGV5m+GcZBFRN5+CiQiCsLXDOMgiyjhfL7GDBFRKT5nGAdZREmnhRNHiYhiyeMM4yCLiLyd94uIKAhfMyw2g6xRTWPv6FvdLVm1m75v0M6ZlTts55bKzrDdhb0529XSJv2mdqltYEGj9JjaarULnl1v5wPLqu26mXmB7bjoGbSdPSN5+7uvDc8ytYOXjp9baspeu22jthkEbdOPm9ofLfq5qa3uP9vUto7MMLW07DW1qxpst0q/2trOUbuvs9KjpvajYxeY2vtbNpravDr7XsmqfTzXXGLNKfueOpK13UPtzbZDa+KfnGud67O21payH9Muqne8aBVQ+Hs+Q62YmF9AuBmWpPwC/MmwpOfXqdbLDHtLbAZZRBQVfztziIhK8zfDOMgiIm/PZyAiCsLXDOMgi4i8PdRORBSErxnGQRZRwqn6G1BERKX4nGGRTvYjInNF5BkR2Sgi60Xk08X6F0Vkt4isKd5ujnI7iOj08ipl3ZKA+UUUH77mV9RHskYB/KmqrhaRVgAviciTxZ/do6pfC/pAU1N53DLl2Mn7D/dPM8s019l5qjrm23mqjvbbbpX1g3NMbX/zm6bWnbZdGBdm7OPVY7tdr9qujve2rTW1RrH7sXvU7u+uETt/12jejpt3NI3/3f4dLWaZlGNKqtFDdl+P5203yO91/cLUvnvwalN7bdDODzZl2kpT2z1qu3p2Zu2+vti3wP5uv+1iun6Rna/NZX9+0NT+Y3Ceqbk6qq6Yss3UljWU/gzj6v55ov9iU3tfyzpTCzp/WRC+ns9QZZHlFxBuhiUpvwB/Mizp+QUww0qJdJClqj0Aeor/PiYiGwHMjnKdRFQ+Xw+1VxPziyg+fM2wSL8uHEtE5gO4HMCJof+nRORVEfm2iNiPKoXfWS4iq0Rk1cFD9tpWRHTmFALV8m5Jw/wi8le5GTaZJmWQJSItAB4G8BlVPQrgmwDOBbAEhU+KX3f9nqquUNWlqrq0q9Me4iSicGiZtyRhfhH5z9f8iry7UEQyKATUv6rqjwBAVfeN+fm3ADwW9XYQ0Sl43JlTbcwvohjwOMOi7i4UAPcC2Kiq/zCmPvbswdsB2DPiiGjy8FCWwfwiihFP8yvqI1lXA/gYgLUisqZY+wsAd4jIEhR2dzuAT5T7wIdytrvkotY9pvbODtuZMb3Ozq315kiXqT09YLs/FjfYdVxeb5/GhRm7fUcd3R9DeszUGsXR/VF/xNSeT9u5yV5oOMfUZswev9ye9jazzNZtM01NBux2pJE3Ndf+/0vezi02mLXzq12YseeqXNFgnxPA1trTdm6txk7b2bQo02Rqecd+uDqv1hy33Tm7h9pN7fK2Haa2sf6AqV2UGf8cuLpr/kvby6Y2NWVfC9c+VCrsT4Ei8icA/isKf+NrAfwegH8G8E4AfcXFPq6qa9yP4IXI8gsIN8OSlF+APxmW9PwCmGGlMizq7sJfAM6psf8tyvUSUXnCbH8WkdkA/hjAYlUdFJGHAHyk+OP/pqo/DG9t0WF+EcWHrxnGK74TJVxEM9jXAWgSkSyAZgD2EAoRUQh8zrBJu4QDEXlKAaiUdwO6TlyeoHhbfvLhVHcD+BqAHSh03/Wp6s+KP/5y8dIH94iIvSokEVG5ys2w0+QXEG6G8UgWEVVyqP2gqi51/aB43ahbASwA0AvgByLyOwA+D2AvgHoAKwB8DsDfVLjJREQnlZlhp8wvINwM45EsIgq7u/DdALap6gFVzQL4EYBfU9UeLRgG8B0AyyLYEyJKonC7C0PLsNgeydqftd0l0+qOm9qNU14zNVeXw73Ddr6p5/sWmtrKlO1+6Z32kqnd0GS7TlpS9sii7eEJzjUf1r4B+4gzm8d351zRucv+3tFWUxvosY+1L2vn1nJ1l3ywc5Wp/c8d7za1D2/+oKl1NNium6vb3zC132nbbGqu59i1fb8atl08VzTY98Xfz7RdMi6vZ23n1QO99u9veccL4+67OoKaxZ5bsGvU1qbaqeQqFPpVkHcAuFJEmgEMArgBwCoR6VbVnuKlEW5Dwi99EGaGJSm/AH8yLOn5BTDDSj1QbAdZRBSiEDtzVHWliPwQwGoUJll+GYVD64+LyHQUOvbWAPhkeGslokTzNMM4yCJKugiulqyqXwDwhQnl60NdCRER4HWGcZBFRCc6boiI4snTDOMgi4gSM1UOEdUoTzOMgywi8jagiIgC8TTDYjvIuqn1VVMbUdtdMb/OdpysHbGdFI/1XGxqDelRU8s7DknuaZ1mapuz201tZtpun6ubJKj1A7NN7fiIbdc4nB4/91Vz3bBZZjRnt01abQdLd32vqbnmn3LND3Zu60FT++mmxfbxRm03zcACu1+XNNoOo2saTcm5fedn7Gu7P2eflxnpYHN1OecXc7xXPv76HePuz2uxc7q9Z9oGU7ukYbephebEhfxoUoWZYUnKL8CfDEt6fgHMsFJiO8giovCEOe8XEdFk8zXDzuhipCLynrA2hIiqKNyLkcYC84uohniaX2d6xfd7Q9kKIqqu8ucurAXML6Ja4Wl+lfy6UEQePdWPANhL9hJR7EiNHJ2aiPlFlAy+ZliQc7KuAfA7ACbObSDg3GNE8VdDXwE6ML+Iap3HGRZkkPUCgAFVfXbiD0RkU/ibdCqKnL7VZeGap8nF1ZkxpHa3bznLdvr0jdrOjNVH5praymN2PrDWlO1gPJqxXRj/3r/I1N7v6Dq6MGPbTrrr+0yto9nOm3V4YPx+1KfsvGQNjm6V7KB9nrYO2TnSMHWnKV3XOGRq26fY5S57m60912ufk/lNh0zN1Ynjsn7Edt1MTdn9Pebo7JqVtt+ou95Tro6droydc+2NLd3j7h+dbXfiAx1rTO2cOrttrnVWpqa+ApzIy/wCws2wJOUX4E+GJT2/AGZYKSUHWap602l+dm24m0NEVeHpp8AzxfwiSghPMyysYSRE5JdhPRYRTbIEdheOxfwiijlP8yvM62QFPPhJRN6pwYFTmZhfRHHmaYaFOcjydBeJ6LQ8vlryJGJ+EcWVxxkW2teFLiIyV0SeEZGNIrJeRD5drHeIyJMisrn4XzuvAxFNGtHybknA/CKKD1/zK8wjWa5h5CiAP1XV1SLSCuAlEXkSwMcBPKWqXxGRuwHcDeBzpR4+LacfE7o6FYbVzl3VnLK1z0zbYmq/tM0l6BmeamrrjtiOi53Hbe7u7bdzkA045up6x2V2Wy7K2HfGH03ban+3+Q1T+9mxS8bdz+Ztl8dQzr4V5rXbbqKc49PCwZztCOpyzJmVEdsR9KGWHaZ219Rgc1zZHhngqUH7fO4ftXOk7RjpstvSttqxDruWvrx9Y+yxu4bFjrnJuuceHnf/7dPt/t8y5ZippVD5HHGBJGTgdBpVzy+g8gxLUn4B/mRY0vMLYIaVUvaRLBFpK36S6xCRjjE/+tjEZVW1R1VXF/99DMBGALMB3Arg/uJi9wO4rewtJyIqE/OLiCZT4CNZIvIJAH8DYBBvjRkVwDkAoKrrSvz+fACXA1gJYKaq9hR/r0dEHBcuAURkOYDlADBvNueyJopKrX8FyPwiqm2+Zlg5f/l/BuAiVT1Y7kpEpAXAwwA+o6pHRYKdoKaqKwCsAIArLmv09CkkqgGenjQaIuYXUS3zNMPK+UhwbesAAB/2SURBVLpwCwD7xXUJIpJBIaD+VVV/VCzvE5Hu4s+7Aewv93GJKCTlXiMrnsMF5hdRrfI4v8o5kvV5AP8pIisBnLzOv6r+8al+QQof+e4FsFFV/2HMjx4FcCeArxT/+0g5G01EIYvnwKkczC+iWuZphpUzyPoXAE8DWAt3Y4TL1SicULpWRE5MZvQXKITTQyJyF4AdAD5cxnYAcHdITE3Z6wlmxHahXJgJNodSZ/q4qfUMtJnajt2208PJ8aXx1Gl2HXtH200tDzv3lcuyBrsfVzZsHHf/UN6uc3u77UwZytu3xyN9V5jaKyN2e29oGjG1323ba2o7Rm3n1HzH6+ji6pxxrfcB2+iCB7e+zdSG5mdMrbu+19QaxW7zs73nm9q2Yx2mtqRrfMfO3856ziyTQpOpufY1TL6ezxAir/ILCDfDkpRfgD8ZlvT8AphhpZQzyBpV1c+W8+Cq+gu4W6MB4IZyHouIIuRpQIWI+UVUyzzNsHLOyXpGRJaLSPcpWqCJKK48PqchJMwvolrmaX6VcyTro8X/fn5M7WQLNBHFU0Ku4s78IqpRPmdY4EGWqi6IckOIqIo8bX8OC/OLqMZ5mmElB1kicr2qPi0iH3T9fExbMxHFlaefAs8U84soITzNsCBHsq5FoSvnAyjshkz4b1VCytWFcyayOmpq59bZuZbqUrZDIr3fdnW4ZM7pN7W3z9ppalc22rmgto3aUfo/HninqV0yxXZ//H7bvnH3O1NTzDKddsosZzfIVY0vmZprjjTXfGtNYp/P2elgXSg/GWgxtaG8fd4/3NJnalc02udYHZ96frD5clPrarOvmcv+PjuvW9rxXrlwds+4+y2piOfzCkiibfypJi/zCwg3w5KUX4A/GZb0/AKYYaUEGWQdE5HPAliHt8IJ8HbcSERl8fh8hhAwv4hqnccZFqS7sAVAK4ArAPw/ALoBnAXgkwAWR7dpRDRpQu4uFJE/EZH1IrJORB4QkUYRWSAiK0Vks4g8KCKOYw+hY34RJUHI3YVhZVjJQZaqfklVvwSgC8DbVPXPVPVPUQitOcE2l4i8FuIgS0RmA/hjAEtV9WIAaQAfAfBVAPeo6kIARwDcFcGejMP8IkqIcD8khpZh5Vwnax6AsZejHQEwv4zfJyJPnWiBDnoLoA5Ak4jUAWgG0APgegA/LP78fgC3RbEvp8D8IqphIecXEFKGlXOdrO8BeFFEfozCWPD24kqIiE5S1d0i8jUUppwZBPAzAC8B6FU9eYb2LgCzJ3GzmF9EFEiYGVbOdbK+LCKPA7imWPo9VX25rC2P2EDezvvUnAp22kda7EE9V3fJwtYDprY247ieoWO4nM/bjpD6lO0K+tWQ/Rbjf+/6dVPb/OYsU2u5fNjUhlvHd+y8nrXrdMk42jUWZWw3zVXOJqmAHUti34Ku7pydI/bi3N/ZdpWpbZy7ztT+qmu9qd109gZT+8FaOx/Yrv5OU2tosc9xdk+zqY1k7HvgyELbFTWRa/9fzw6amuu1qFj5J412iciqMfdXqOoKABCRaQBuBbAAQC+AHwC4KZS1VigO+QVUnmFJyi/AnwxjfrnFIMNOmV9AuBlWzpEsqOpqAHYWTiKKr8o6cw6q6tJT/OzdALap6gEAEJEfAfg1AO0iUlf8JDgHwJ4Kt7gizC+iGlV+hp0uv4AQM6ycc7KIqFaF2124A8CVItIsIoLCZMobADwD4EPFZe4E8Eio+0BEyRVud2FoGcZBFhGFOshS1ZUonBy6GsBaFHJmBYDPAfisiLwBoBPAveHvCBElUoiDrDAzrKyvC4mo9gjCv5Cfqn4BwBcmlLcCWBbumogo6XzOMA6yiIjXPyeiePM0w2IzyBIAqTHfbrq6FxocnR45dczV5ejEcS339OBUU0s51rtwiZ2rq3/Ezuc0o9nOI7XhiO2weX7PAvt4m6aZGjpt99CRUdsl8qHNHxh3/zxHh9Hvd/7C1FydH7tGj5va9LTtfurL2207kEub2kX19nlKOb7Ffl/La6b2rdzVpvb8QUenlKM7Z9vxLlOTA3ZbGo46ZnZXu79tR+xigzMdm1J3bNz9h/vt6/qLo4tMbVnrFlO7IHPYrqASlZ34TmWYmF9AuBmWpPwC/MmwpOcXwAwrJTaDLCKKkKcBRUQUiKcZxkEWEXkbUEREgXiaYRxkEZG3h9qJiILwNcM4yCIibz8FEhEF4mmGRXqdLBH5tojsF5F1Y2pfFJHdIrKmeLs5ym0gohLKvUaWp2EWBWYYUQx4nF9RH8m6D8A/AvjuhPo9qvq1iNddliP5IVO7osHWrmw8ZGor6i43tQ393ab2m12rTC2r9iX40rrfMLX0sO0SyZkK8MohO19lOjW+o+htLW+aZfbmWk3thb7pprbLMQfXu1rtPFp7R23X0U8OXWpqrrnPPtr1gl1Hk+06+rPznzS1xw4uMbUPvG6fz9d2OuZN222f44xtpkHDUduhdWyu/bwy0mlfoed7zxt3v3fEdj81pu1z8lsdB+2GhPgZyddD7R64DzHIsCTlF+BPhiU9vwBmWCmRHslS1ecAhNSjSUSR8fiTYDUxw4hiwtP8qta0Op8SkVeLh+IdF1ApEJHlIrJKRFYdOOT63ENEYRAt70alM4z5RTR5fM2vagyyvgngXABLAPQA+PqpFlTVFaq6VFWXTu+0F7IkopDwSFY5AmUY84toEnmaX5M+yFLVfaqaU9U8gG+Bc5kRVRdPfC8LM4zIMx7n16QPskRk7BmVtwNYd6pliSh6UsEtyZhhRH7xOb8i7S4UkQcAXAegS0R2oTCj9XUisgSF8eR2AJ8I8lgK91xfpYy6elgcI9lpqcayH/uE5tSIqc1vsl08U1LDpnZDk23/OHDBc6a2aZ7t9tnQZyeXSjm+cB7Njx9Lf/nVm8wydWn7PM1otXOVzXTMX9abtZ0z63pt98uuQ/bUlaZG+9xd177J1PrzO+1yzbbDaEurfU6+u/YdpobDdv6unOMtMOw42+Zoxn420ZR93hv226+IVj954bj7Q7Pt/Gifu/pxU7uq8rdnMAk/OnUqYWVYpfkFBMuwJOUX4E+GJT2/AGZYKZEOslT1Dkf53ijXSUTl48nsbswwonjwNcN4xXci8vZTIBFRIJ5mGAdZRORtQBERBeJphnGQRZR0CkhlpwsREVWfxxnGQRYReXs+AxFREL5mWE0NstLiuCKF44l/cdgu5+p8yKkdGrvWMTPTF6i2tOGoXQnsij/ZvsPU+tpeN7VV09pMbcPQHFN7ofeccfe3bbFdM9khu1/bmm3HzYGuFlNra2o3tZ6tXaaGBvt83r7oFVNrTw+Y2r6c7RxaP2L3Y0bGPsfdXfa16HvZ/q6rO0cXHTe1c2fYObhe23KWqeWP2+6cTN/451ky9jl5eI+dS277kH0+396y1dQq5mlAJU2lGZak/AL8ybCk5xfADCulpgZZRFQZXz8FEhEF4WuGcZBFlHS8ijsRxZnHGcZBFhF5G1BERIF4mmEcZBElnMDfQ+1ERKX4nGEcZBGRt58CiYgC8TTDEjnIWubocnPNlZ0R+/S45h9rT9sOjs6UrU0NOL9YyrEtrt+9ocnOm3VZvZ2rtnnCnGOjF9nH7x22j//GNtvB0n9giqk1nDVqak6Ors7H31xsaq9MnW1qS6fZjqXtg52mtneg1dQO9NluogbHdHCOKdwwr+uwqX1q7tOm9n+n2G6aQyP2uTo81DTuft+E+wCwZdcMUxvI2rnKzm3cb2qVEvU0ocjJZlhy8gvwJ8OSnl8AM6yURA6yiGgMj08aJSIqyeMM4yCLiLw9n4GIKAhfM4yDLCLy9lMgEVEgnmYYB1lEFPqnQBE5H8CDY0rnAPhrAO0A/gDAgWL9L1T138JdOxElTZgZFmZ+cZBFRKF/ClTVTQCWAICIpAHsBvBjAL8H4B5V/Vq4aySiRAsxw8LMr9gMskaRx8HcW3NCdaXtnFRBuebvcnXEuDpxjuQHTe3XG2xbR3PKdlIE5VpvUK7n5ZLGnePuz5pl58I66pj46htDN5jaoa0dpta/xnbJNNqGIFcDFI4PTjW1N7obTG3bAbsOl5H9ttOlaa+dg6uh1/5FZqeIqR0dss/LBfV27q8rup80tfUjdt/SE6aKXzlwrllm24zppjbTMW/c/PoDplYRjfx8hhsAbFHVN0Xsc5wEE/MLCDfDkpRfgD8ZlvT8AhKRYWeUX47/7RFR4miZt/J8BMADY+5/SkReFZFvi8i0M9twIiJ4m18cZBEl3ImrJZdzA9AlIqvG3JY7H1ukHsAtAH5QLH0TwLkoHIrvAfD1yHeQiGpauRmGScyv2HxdSEQRKv9CfgdVdWmA5W4CsFpV9xVWU/gvAIjItwA8Vu6KiYiM8jJs0vKLR7KIqJIjWUHdgTGH2kWke8zPbgdgL/FNRFQmX/Mr0kFW8TvL/SKybkytQ0SeFJHNxf/ynAyiair3fKyAISUizQDeA+BHY8r/XUTWisirAN4F4E9C2YeIMMOIYsDj/Ir668L7APwjgO+Oqd0N4ClV/YqI3F28/7lSD1QHwbQx81/ltPIOFtfvphwdhy7NjvnA6mC7P3yyrGFiJ9KAWeaFYdv50XvUdvqkh213RSpr15l3NCdl+m1NO+07vskxp9nxo7ZLpnWq7ZQaabKvrabt6zM81e5Hzjb2YHDAFtcOzzS1W6YcM7WMo5umTcZ3Hl3d8JpZZkfbGlMbUvv+nJ6q/G9gIkfT0BlT1QEAnRNqHwt/TZG6DyFk2MT8AsLNsCTlF+BPhiU9v4DazbCw8ivSI1mq+hyAiTNU3grg/uK/7wdwW5TbQEQBRNtdGFvMMKKY8DS/qnHi+0xV7QEAVe0RETtdNxFNKl/n/fIUM4zIM75mmNfdhcW2yuUAMHe234e0iWJLUUl3IZXA/CKaJB5nWDW6C/edOEO/+N/9p1pQVVeo6lJVXTq9kyFFFJUIuwtrUaAMY34RTR5f86sag6xHAdxZ/PedAB6pwjYQ0Vg8J6sczDAi33iaX5F+XSgiDwC4DoWrq+4C8AUAXwHwkIjcBWAHgA8HfDTnnINjubp1XrKNHtgyYk+huKPVzufkXod9herOYDq2vvyQqU1N2U6UifOeAUDG8XxMS9lumolzibnmFtuZtXNrvW/hRlPb1m2Xa66zrTlbjtjlBtbYOcNcb/hzOg6Z2vy5E889Bs5r3mdqP9u/2NTWp2abmn3W4fyIM6vFPu8HRttM7WDObsuMdIupBZnXbUFd5XPaVeLE1ZLJCi/DSucXUHmGJSm/TlWrSoYxv5yYYW+JdJClqnec4kd21k4iqg5Vb89nqDZmGFEMeJxhXp/4TkSTI4rrZBERTRZfM4yDLCLy9lA7EVEQvmYYB1lESacA8p4mFBFRKR5nGAdZRMSOQSKKN08zrKYGWa7unbl1x01t80jlV65ocMz9tTtn56CanbZzRrm2z9WJ49Kayji2xdaCcHX6LGmwtd8860hFjw8AH97yPlN7c5ftLsxOsc/JzKvsPFr3nLXS1B7ut/Py7j421dRmzOwztQs7bDdNyvGl/sv755jaf/adZ2rXNL9h1xujSyP5eqg9aaLOsFrIL8CfDGN++cPXDKupQRYRVcjTzhwiokA8zTAOsojI20+BRERB+JphHGQRJR2v4k5EceZxhnGQRZRwhasle5pQREQl+JxhHGQREQLOlkFE5CdPM6zmB1nd6Smm9vbGNx1L2m4alx05O2tUR8q2Ybg6dubV2W1JOebods0PVYdgrR6u3x3W8fNyueYMOzdju4Rc2+ayO2e7aV4/ON3Uph527NegnTht94DtsHG5oH6vqX14/st2vWnbdbSwwf7uI0feZmrdrUdNbWaDrQ2pfX2O5O16+/K5cffn17WaZYLODxb09QnC10+BFDTDkpNfgD8ZlvT8AphhpdT8IIuISvD4fAYiopI8zjAOsogSz9/JVYmISvM3wzjIIiJv25+JiILwNcM4yCIibz8FEhEF4mmGcZBFlHQKOGbkICKKB48zLJGDrEWZYJ04ObWv2nRHJ06TYw6uFkczjevxYJtT3B0XjuVcj+eaX2yioPONBe0Q2Zq13SXHd7eYWn2n3ba6QfvpY/Ne29XTf96wqV1U32Brna+dcjtLefKo7WI6nq03tYzYDpvW1Kipbcra57k33zyhcsgscyBnf68jbfe/Mczj455+CiS3IBlWq/kF+JNhzC+AGXZ6iRxkEdEEfuYTEVEwnmYYB1lE5O01ZoiIgvA1wzjIIiJvD7UTEQXiaYZxkEWUdApvp6QgIirJ4wyr2iBLRLYDOAYgB2BUVZdWa1uIkkyg3h5q9xkzjMgPPmdYtY9kvUtVD4b1YFuydq6l5pR94menK59rqTdvuzBa6myXiIurmyZox45L0E6cBkf3UKV+aac+w+NHLzW12QsPmFr/PNvpknO8Pu+dtd3UVvRebGoLGvab2u1T+uwGOuwaPW5qFzTtMbW5s23nzEDe7scXd/+Gqb3cM9vUcvnxr9nIgH1tGpptl9C0FvvePm/qqf507j1F/TRCDigROR/Ag2NK5wD4awDfLdbnA9gO4LdU9UioK59cXmQY8yu4MDMs6fkFRJFhFeQXEGqGhZlf4c3OSETxpVrereTD6SZVXaKqSwBcAWAAwI8B3A3gKVVdCOCp4n0iojPjaX5Vc5ClAH4mIi+JyPIqbgdRsp04n6GcW3luALBFVd8EcCuA+4v1+wHcdmYbX1XMMCIflJth5Tmj/Krm14VXq+oeEZkB4EkReU1Vnxu7QDG4lgPAvNnV/maTqHZFfD7DRwA8UPz3TFXtAQBV7Sn+/cfVaTOM+UU0eSLMsDPKr6odyVLVPcX/7kfhMNwyxzIrVHWpqi7t6nRcgpiIwlH+14VdIrJqzM15JEdE6gHcAuAHk7k7k6FUhjG/iCaRp/lVlY9XIjIFQEpVjxX//V4Af1ONbSGiYOcpTHAwYDfdTQBWq+q+4v19ItJd/BTYDcCe/RsDzDAin5SdYZOWX9U6hj0TwI9F5MQ2/H+q+tMzfdDjanfnuJ2mCd3pYF/KurpfNmQ7TW1OnaNdJaCe3KCp7XPM+3Rhxu5Ic8p2ibhMnEvM1YkUdB6xZQ12uefSdv+Xz/8PU5ufsZ0kHWnbcTI9ZdexJ2df2zVD80xtf67H1Gal7RxkU1J2365vfsPUpjrment20M5N9i8vXWtqHc/brq3M8fFBoI4DHP1n2dd/b/cUU/ujm56xvwzg/zirp6EAcpEdar8Dbx1qB4BHAdwJ4CvF/z4S1Yoj5n2G1Wp+Af5kWNLzCwg/w8rOLyDKDDvj/KrKIEtVtwK4rBrrJiIrivMZRKQZwHsAfGJM+SsAHhKRuwDsAPDh0Fc8CZhhRH4JO8PCyi+ejUlEkUxJoaoDADon1A6h0K1DRBSekDMsrPziIIso6RRA3s+rJRMRleRxhnGQRZR4FZ34TkTkCX8zjIMsIvI2oIiIAvE0w2I7yOrL226QRrFPcodjbqmco4Mn6Dxa722y63XO3xVQe8puy2P9C0xtenqjqc1zdOdsG7XdLhOfl+607fJw7X9/ftjU/m3gLFPrzvSa2m+3BuvMP5izz92A44/l0nq7r7PSW00tI3bitKza+dqmpZocNbt9Lw7b7Tucc3T7tNn3Rf0xO6fXtOd3jbs/ssB1LTvb1TMw327Hzc12rrKKeRpQtSzMDEtSfgH+ZFjS8wtghpUS20EWEYXE4/MZiIhK8jjDOMgiSjwFzuBoBhFRdfmbYRxkEZG3h9qJiALxNMM4yCJKOo8PtRMRleRxhnGQRUTefgokIgrE0wyL0SBLx3XBpGG7MBrFzo/V5ehCcdmYDTZ/10WZZlu0m4L9Odsl05my8zm1pGwXxu+22TmoMmJ/12WKY1s6HZ0oE7k6jBrEvj1+c8oBU3th2C6XVftaNIjtVulK2+dz/YjtCHJxTKOFYzk7kdaS+mBvc9d8aP/n0K+Z2hOvLza10X67b62OTqHRHeO7czJt9v3Z1Gz3ITVg92HA8RxXzNOAqh1q/sbCzLCk5xdQnQxLen4BzLBSYjTIIqJo+HshPyKi0vzNMA6yiJJOAeT97MwhIirJ4wzjIIuIvP0USEQUiKcZxkEWEXkbUEREgXiaYRxkESWeetv+TERUmr8ZFqNBloybm6pFbFdLk6N7LahFdXZuKddcWK55pFzLzXB0nQSdI6zZMadXUEFmMAu6HRmxHSIuZ9cdM7Vhx/yQQR9vejpYx8mCOvscv54dDPS7rtdxFHa9adg/3PxB+95r226f+Xyd/d3hDywbd19TtoOnd4F9nlrOPWJqW7KOTrFKKKCeXi25dojJiTAzLEn5Vc62VCPDkpRfADOslBgNsogoMp5+CiQiCsTTDOMgi4i8PZ+BiCgQTzOMgyyipFP1tv2ZiKgkjzOMgywi8vZTIBFRIJ5mWNDzDEMnIjeKyCYReUNE7q7WdhARoPl8WTdihhH5xNf8qsqRLBFJA/h/AbwHwC4AvxKRR1V1w5k8rqtLxtWF4lrOVQu6DpeUa/zqmJcraJeMa72udbjm0jJzpgXch6Dm1Nm5qw465j5rCvhauDqbXFz7f0HGbotrTq9/6j3P1FrTtrNnMG+7vbTFdvY0HbTLjTbZF3z3jdnxj7XPMafbdLsdf3fhE6Z2SX2wOdJK83dKCl/5lmFJyq9TPd6ZqDTDkp5fADOslGodyVoG4A1V3aqqIwC+D+DWKm0LUbIpCp055dyIGUbki3IzbBJV65ys2QB2jrm/C8A7qrQtROTpNWY8xgwj8omnGVatQZbjwLO9apqILAewHADmzeY5+kRRUFVoLtgFFOmkkhnG/CKaHD5nWLW+LtwFYO6Y+3MA7Jm4kKquUNWlqrq0qzPYlXuJqHya17JuVDrDmF9Ek8fX/KrWx6tfAVgoIgsA7AbwEQAfrdK2EJGnh9o9xgwj8omnGSZapTPyReRmAN8AkAbwbVX9conlDwB4E0AXgIPRb2HkamE/amEfgNrYj7H7cLaqTg/6iyLy0+Lvl+Ogqt5Y5u/UlHIybEx+AbX3fouzWtiPWtgH4K39KCu/gIoybNLyq2qDrEqJyCpVXVrt7ThTtbAftbAPQG3sRy3sQ1LUwmtVC/sA1MZ+1MI+ALWzHxNV7WKkRERERLWMgywiIiKiCMRxkLWi2hsQklrYj1rYB6A29qMW9iEpauG1qoV9AGpjP2phH4Da2Y9xYndOFhEREVEcxPFIFhEREZH3YjXIiuus9yLybRHZLyLrxtQ6RORJEdlc/O+0am5jKSIyV0SeEZGNIrJeRD5drMdmP0SkUUReFJFXivvwpWJ9gYisLO7DgyJSX+1tLUVE0iLysog8Vrwfu31IGuZX9dRCfgHMsDiKzSBrzKz3NwFYDOAOEVlc3a0K7D4AE6/JcTeAp1R1IYCnivd9NgrgT1X1QgBXAvjD4vMfp/0YBnC9ql4GYAmAG0XkSgBfBXBPcR+OALiritsY1KcBbBxzP477kBjMr6qrhfwCmGGxE5tBFmI8672qPgfg8ITyrQDuL/77fgC3TepGlUlVe1R1dfHfx1D445iNGO2HFvQX72aKNwVwPYAfFute7wMAiMgcAL8B4H8X7wtitg8JxPyqolrIL4AZFkdxGmS5Zr2fXaVtCcNMVe0BCgEAYEaVtycwEZkP4HIAKxGz/Sgeol4DYD+AJwFsAdCrqqPFReLwvvoGgD8HcGIeiU7Ebx+ShvnliTjnF8AMi5s4DbJKznpP0RORFgAPA/iMqh6t9vaUS1VzqroEhQl9lwG40LXY5G5VcCLyfgD7VfWlsWXHot7uQ0LxNfJA3PMLYIbFTbUmiK5EyVnvY2afiHSrao+IdKPwqcRrIpJBIaD+VVV/VCzHbj8AQFV7ReTnKJyf0S4idcVPUb6/r64GcEtx3rxGAG0ofCqM0z4kEfOrymopvwBmWFzE6UjWyVnvi10HHwHwaJW36Uw8CuDO4r/vBPBIFbelpOJ35vcC2Kiq/zDmR7HZDxGZLiLtxX83AXg3CudmPAPgQ8XFvN4HVf28qs5R1fko/A08raq/jRjtQ0Ixv6qoFvILYIbFkqrG5gbgZgCvo/Ad9F9We3vK2O4HAPQAyKLwifYuFL6DfgrA5uJ/O6q9nSX24ddROHz7KoA1xdvNcdoPAJcCeLm4D+sA/HWxfg6AFwG8AeAHABqqva0B9+c6AI/FeR+SdGN+VXUfYp9fxf1ghsXsxiu+ExEREUUgTl8XEhEREcUGB1lEREREEeAgi4iIiCgCHGQRERERRYCDLCIiIqIIcJBFREREFIE4XfGdJomIfBFAPwpX4n1OVf/9DB5rCYBvFh8rB+DLqvpgGNtJRDQR84t8wkEWnZKq/nUIDzMA4HdVdbOInAXgJRF5QlV7Q3hsIiIn5hf5gF8XEgBARP5SRDaJyL8DOL9Yu09EPlT893YR+XsR+aWIrBKRt4nIEyKyRUQ+earHVdXXVXVz8d97UJgbbPok7BIRJQTzi3zFI1kEEbkChTmkLkfhPbEawEuORXeq6lUicg+A+1CY6LMRwHoA/xxgPcsA1KMwrQgR0RljfpHPOMgiALgGwI9VdQAARORUE9eeqK8F0KKqxwAcE5EhEWk/3SH04gz33wNwp6rmQ9x2Iko25hd5i18X0glBJrEcLv43P+bfJ+6fcsAuIm0AfgLgr1T1hYq3kIjIjflFXuIgiwDgOQC3i0iTiLQC+EBYDywi9QB+DOC7qvqDsB6XiKiI+UXe4teFBFVdLSIPAlgD4E0A/xHiw/8WgGsBdIrIx4u1j6vqmhDXQUQJxfwin4lqkKOsRERERFQOfl1IREREFAF+XUihEJFLUOi+GWtYVd9Rje0hIgqK+UVR4deFRERERBHg14VEREREEeAgi4iIiCgCHGQRERERRYCDLCIiIqIIcJBFREREFIH/H6PNfUOHNnjcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=[10, 10])\n",
    "ax1 = plt.subplot(221)\n",
    "ax2 = plt.subplot(222)\n",
    "xds1[0,:,:].plot(ax=ax1)\n",
    "xds2[0,:,:].plot(ax=ax2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f2eae189050>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAEmCAYAAABVkVXzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5gcVZk/8O87k8n9ngkwJIGEEC6KEnCIAZRrYAFB0BUXXV1Q1shvZQXF3wKLi7Aru7ir4AVWf1FY4j6sKLcH5CLXIItyCyFAYoAkQMiQkBCSyYUkc31/f3QHZ+Y9PVPVXdV1Ttf38zz9kHnTXXVqZvLldHW9dURVQURERETpqMt6AERERES1jJMtIiIiohRxskVERESUIk62iIiIiFLEyRYRERFRijjZIiIiIkpRVSZbIlIvIs+LyD3Fr6eJyNMislxEfi0ig6sxDiKiuJhfRFQpqcZ9tkTkmwCaAYxW1VNF5DcA7lDVW0TkZwBeUNWf9reNIaPG6sjGpve/njpuWKJjbG3rMrVBdWJqIxuSnZ++19FtaivefMc+Ue3zDt63yT4vgp1d9me+s9Nu3x69W2e33d6EYYMivXb15jZT27xlp6nV1dvv+4F7jjI115jtkbmfF/Vfgus34IXX1pta/SD7/+CD9ho74PZd4920o9PU1rfucL5+x9srNqjqxAF3VDRFhulO515L24D2B1T1pFgvClQa+QUkm2F5yi/AnwxjfrlVkmFx8wuIn2FZ5Fe0/yNWQEQmA/gEgKsAfFNEBMBxAD5ffMp8AFcA6DesRjY24S+u/OX7X9945kGJjvPeFZtMbdywBlM7YtLIRPf77NrtpnbG+T8zta52+0v5h99eVtY+V7Z2mNrSd7aZWoMjrF3e3d5ual/40G6RXvuNe181tft/97KpjRpvv++/v+I4U2twDLnN8W/QdWiOvHUa4kiryZ/9kamNmLiXqf3h+k8NuH3X/0hu+5P9H9hP7ljqfP2L3/vEqgF30kMbuvFZifc/vv/UVY2xXhCotPILSDbD8pRfgD8ZxvxyqyTD4uYXED/Dssiv1CdbAH4I4B8A7JrGTwDQqqq7prktACZVYRxE5CAA6iXqOYCi/Cw8wfwi8lzsDMsgv1K9ZktETgWwXlWf61l2PNV56CIyV0QWisjCnVtbUxkjEQH1Eu+RB8wvonD4nl9pn9k6EsAnReQUAEMBjEbhneJYERlUfHc4GcAa14tVdR6AeQAwYdqB+XkvTVRFZZ3ZygfmF1EAQsiwVM9sqeqlqjpZVacCOAvAo6r61wAWAPhM8WlnA7grzXEQUT9intXKy5kt5hdRIALIr2pcs+VyMYBbROS7AJ4HcMNAL1AAXT2uAPzyrUsi7ejMQydHel6dY1Yc9WLS+1cm+xHBL/79y6bWnWDX6MpN9qLWqBeSujz6su1keeYNe8FuV8QrOI+fs5+pnXrQHqbmupjU5Qs3Lza11Y7xPfVPx0bboMPWNStNbcemdY5nDnyB6VBHGrgu1nVcgwoAOO97A+6ilxDeFXqm4vwCks2wPOUX4HeG5T2/gMoyLG5+AWFkWNUmW6r6GIDHin9+DcCsau2biEorBFXWo/Ab84vIXyFkWFZntojIG+L9u0IiotL8zzBOtohyTsB1u4goXCFkGCdbROT9u0Iiov74nmGcbBHlnOSow5CIak8IGRbMZEsA1PfoNona2RZVfcROlqQ7d6IaM9QuveHiWtphc59100YOrrfP2WnXrYrK1S3l+n7euqjF1AYPsmNxvfbk6dHW5Dr4Hx4wtaEj7Bpfj0fs3PnJM2+Z2nXznzE17bZr03V32CVAkvTXJZYTOa+Mbfn+rjB0ffMLSDbD8pRfgD8ZxvyqjCvDyskvwP8MC2ayRUTpCKGTh4iolBAyjJMtopwL4R41RESlhJBhnGwRkffvComI+uN7hnGyRZRzhYtLPU8qIqISQsgwTraIyPt3hURE/fE9wzjZisnVaZJ0Z6SLq/vGZYjjzm7Prdma2Dhcx3/itDGRXhv1eVF95falprbzPbtm2h5TbSeQ6/vkcvMDr5qaq3PHxfW8DTt61xqHRfu5ugxKKFxCuLiUkpH3/AL8ybBayC+AGRYVJ1tEOScpLHUhIkMBPA5gCAo5c5uqfifRnRARIfkMSyO/ONkiojTeFbYBOE5Vt4lIA4AnROR+VX0q8T0RUe4lnGGJ5xcnW0Q5l8bFpaqqALYVv2woPtL/vIqIcifpDEsjv3xfu5GIUrbreoc4DwCNIrKwx2Ou2a5IvYgsBrAewEOq+nR1j4yI8iBuhiGD/OKZLaKcEwCDIy730sMGVW3u7wmq2gVgpoiMBXCniBykqkvKHCYRkVMZGVb1/ApmsqVIt2umkm6Uaqw3NnO3YanvI4qPNI0q+7XqOM37wjrbfdPWaRdIu2vJ26b25B/ftPtwdNAc/+GmqEO0+73oKFM7+uJ7TW2Pg+1aZVLnWL+tz+JvlXTyJCnNe9SoaquIPAbgJAC5nGwxv/zILyDZDMt7fgG1n2FJ5Rc/RiTKORGgTiTWY+BtysTiO0KIyDAAcwC8nPKhEFEOxc2wgbeXfH4Fc2aLiNIikOTbEZsAzBeRehTe1P1GVe9JeidERClkWOL5xckWUd4JUJfwZEtVXwRwSKIbJSJySTjD0sgvTraIck4ASD2vKCCiMIWQYalOtkrdhVVEbgJwNIDNxaeeo6qL0xwLEZUgSONjxOAxv4gCEUCGpX1my3kX1uLf/V9VvS3qhsYPH4yzmqckMqihg6LNgBc5Ok027uhMZAz92X3k4EjPW7x+h6m5un4a+sz4O7psR0lUS9a/Z2odji4r1z76jqPU8+5xdO64zJptfx/+8NhyU/vGkXtF2t6Vj75uamc3Tza1kWNHmFrDkNGm9uNzZ5naXmMaen3d6WhQ63LUHn5tk2NbCXV4iST+MWKN8DK/gGgZVqv5BfiTYXnPL4AZFlWqky3eRZooDFLn9yn4LDC/iMLhe4alPrp+7sJ6lYi8KCLXisiQtMdBRG5SvLg0ziMvmF9E/oubYVlIfbKlql2qOhPAZACzROQgAJcCOADAYQDGA7jY9VoRmbvrdvpbNr2b9lCJckvqJdYjL5hfRGHwPb+qdt5NVVsBPAbgJFVdqwVtAP4LgP1wuPCaeararKrNo8dNqNZQifJFBFJfF+uRN8wvIo/FzLAspLrXUndhFZGmYk0AnIGcLuFB5AMBP0Z0YX4RhSFuhmUh7W5E511YReRREZmIwvdoMYDzUh5HL/uOH25qC1ZtMbWdjjWukubqtDl0dzs+l8mjonX9zJnau9Pkwdc3m+e41m3bd4Idx4p3bYdTVFE7iE49aI9Iz7tl4WpTm3PigaZ20T2vmNpgx/e93TG+a37/mqmJY8HTJ799TKlh9rJqa0evr8/5xbPmOW+/Ybt2ujvbTc21jlpZxH1M5Gd+ATbD8pRfgN8Zlqf8AphhUaXdjei8C6uqHpfmfokoDkFdDj8aHAjziygU/mcY7yBPlHcB3BCQiKikADKMky2inJMAgoqIqJQQMoyTLSLy/hQ8EVF/fM8wTraI8k7yde8sIqoxAWRYLidbXWq7VqrRueMSdZ1Gl7XbbGfI4nV23a8xQ3r/mF1dOyEaOdSu01Xv6EhxHe8ORxfMdWfYTqDLHlhhal867QBTe3ljm6kdMN7eWHzvUb3HPHO/RvOcBeu2mVpHmz2uwcPs8ZdDANR53slDvfXNsDzlF1AbGVYL+QUww6LK5WSLiHoQ5PJGpURUIwLIME62iCg3Nyolotrke4ZxskWUdwFc70BEVFIAGcbJFlHOSQCn4ImISgkhwzjZIiLvT8ETEfXH9wzjZIso7wJYV4yIqKQAMiyXk61l79j24qwMqeDU58YdtnXatVDq9o6BF/t0LdhaiZOnj010ew+/YRfa3dHeaWp7Txhhapcfv4+pXXD3y6Z257INpvavJ043tcXrd5jauT97ytQeufgoUxva593X1SfvZ54z5BO2trPLtn+LuMNl7A+d5ZIkgHXFqDdfMsyX/AL8zrBazS8g+QyLm19AGBmWy8kWEfVQB9QNrs96FERE5QkgwzjZIso9gdT5/a6QiKg0/zOMky2inBNJfl0xEZkC4JcA9gDQDWCeqv4o0Z0QESH5DEsjvzjZIso9SaNtuhPARaq6SERGAXhORB5S1T8lvSMiyrvEMyzx/OJkiyjvUrhHjaquBbC2+OetIrIMwCQAnGwRUbISzrA08iuYyVZ9nWDM0D8P17UOaYOj9dPV8VKJ4Q3RLsKL2kGzz7ihkZ63brvd3p6j7UKhvnjw9c2m5lpQ1bXw6onTxpjanKmjHXvZy1TufmmtqYlj4fEfn7a/qV38O7to66cPmGBqH97Ndj3desGRpjbY+Y+/91iGRMwHVxdQctK93kFEpgI4BMDTqe3Ec33zC8gmw5hf0UXJsLznF1D7GZZUfgUz2SKilAgg9bE7eRpFZGGPr+ep6jyzaZGRAG4HcKGq2v53IqJKxc+wqucXJ1tEOSflXe+wQVWb+92uSAMKQXWzqt5R7viIiPpTRoZVPb842SLKOwHqEj4FL4W7Fd4AYJmqXpPoxomIeko4w9LIr1RvTCEiQ0XkGRF5QUSWisiVxfo0EXlaRJaLyK9FZHCa4yCi/kl9XaxHBEcC+CKA40RkcfFxSrpHkSzmF1E4fM+vtM9stQE4TlW3FU/JPSEi9wP4JoBrVfUWEfkZgHMB/DTlsRCRiyR/6wdVfQKA34uVDYz5RRSChDMsjfxKdbKlqgpgW/HLhuJDARwH4PPF+nwAV2CAsBrZUIcjJo3sd38vvbPT1KJ28kTt0vn6dX+M9Lwfn3+Eqbk6fJast2ucTXV0+LQ71pXymavzsMHxj8G1FpqrCyhqh4+r9o17XzW1ax1rd33vpH1NzaXO0ZGz5wjXP6WBf2bdEf89dzs6kupKrI0YlwDe3305C9XOL6D8DGN+JS9KhuU9vwBmWFSpX7MlIvUAngOwL4DrAawE0Kqqu1bhbEHh/hVElIUUzmzVCuYXUQACyLDUJ1uq2gVgpoiMBXAngANdT3O9VkTmApgLAFOmTEltjES5lsJNTWsF84soAAFkWNVGp6qtAB4DMBvAWBHZNdGbDGBNidfMU9VmVW2e2NhYnYES5VBdfV2sR94wv4j85nt+pd2NOLH4jhAiMgzAHADLACwA8Jni084GcFea4yCi0kQKd1+O88gD5hdRGOJmWBbS/hixCcD84nUPdQB+o6r3iMifANwiIt8F8DwK97Mgooz4fgo+I8wvokD4nmFpdyO+iMKaQn3rrwGYVdG2HV0Ma7e1lb09V5fO8ArW7ora9eNy9z8dV/Zr07bvBLuu1op3t0d67dBB9h/D7iPsLYpattiOrPtXtpraydPHRtrv5u3tpvbDp1pM7cLZkyNtLyrX72jfdc5eb7Vjm79wtam9+559nmtdtrIEcHFpFtLMLyDZDGN+RZdkhuU9vwBmWFS8gzwR5eajQSKqTb5nGCdbRDknIqiLvxA1EZEXQsgwTraIyPtT8ERE/fE9wzjZIsq7AO5RQ0RUUgAZxskWUe6J99c7EBGV5n+GBTvZcnVFuNaycjn3Wz83td32P9jULv0rW/vaP91sauP3mmFq7yx/MdJYXD75L7bm6vDJYr2xqF07Y4baX60oa8MBwIcm2rXV/nOhvW9kp9punkGO5pZhg+1Yzv9o+p07jl9Rs4rYtLG2m2nS+GGmtrWt09Si/r4PRAJ4V1iLksww5ld0aWdYnvILYIZFFexki4gSEkDbNBFRSQFkGCdbROT9KXgiov74nmGcbBHlnQikzu+2aSKikgLIME62iAjwPKiIiPrleYZxskWUdyKQhoasR0FEVJ4AMiyYydaW9i4sWLXl/a87Eupi2GXH1vdMbbJjbbEPHHWYqa1Y+Kqp1Q+x3RjdHXZtKO3uMjVXJ9DhfxOtO+j3/3VBpOdFEbVrx2V7R3di4wCAv2ves+zXXn/6AQmOBFj7nu2q+e7DK0xtsGMtyGs/sV+vr+tgf4/PPaTJ1C7fYH8Wwwcn9U5OvH9XGLq++QUkm2HMLzdfMqxW8wtghkUVzGSLiFIi8D6oiIhKCiDDONkiyjkJ4IaARESlhJBhnGwR5V0A7wqJiEoKIMM42SLKPf+vdyAiKs3/DONki4i8PwVPRNQf3zMsmMmWonf3jmtNpW7XYk4OE6Z/yNTe2/C2qX2g0Xbk3Dl3lqmdP3GEqS347dOmNmTiXqa2dY3tAonK1Ql09Jd+VNa2bvj+V8oeh0tHl+3kefiNLabWvKdda2z1Ztv1dM+ydaY2e+9xpnbbYrsG2fH7TzS1Tx8wwdQ6IzaH7T4iWovxqCED//Ny7fM//neVqU0Yadcf+8YR9vcJAK4aeGi9if/vCkPXN7+AZDMs7/kFZJNhec+vUvutJMNi5xcQRIb5PRUkoiooBlWcx0BbFLlRRNaLyJIqHAAR5VrMDBtoaynkFydbRHkngNTXx3pEcBOAk9IdOBERYmdYBDch4fwK5mNEIkqLAAlf76Cqj4vI1EQ3SkTklGyGpZFfnGwR5V0AbdNERCUFkGGpfowoIlNEZIGILBORpSJyQbF+hYi8JSKLi49T0hwHEfVHIHX1sR4AGkVkYY/H3KyPImnML6JQxMswZJBfaZ/Z6gRwkaouEpFRAJ4TkYeKf3etqn4/6oZGDx6EOdPGvv/1Ays3lT2od1e+ZGrimBX/zrGPM/Ybb2q/ONN2N357jO0EenrZejuWEUNNbYRjTbP1b9huo51b3jG1rrYdppa24Q32e7e9w3Yaubp7rnn8DVNbv7XN1OrrxNRu22SPddvODlN7szXa9+TRNzab2q2LWiK9duRQ2+Fz5QnT7RP7dJsNsoeFFsdxTR5nf59ca5KVLf4p+A2q2pzcALyUWn4ByWYY86sy5WZY3vMLCDbDqp5fqU62VHUtgLXFP28VkWUAJqW5TyKKScT5P+u8Y34RBSKADKtaN2LxYrNDAOy6gcv5IvJiscXS3nCk8Jq5u07zbdhg3wURURJSufXDrwA8CWB/EWkRkXNTP4wUMb+IfJb4rR8Sz6+qTLZEZCSA2wFcqKpbAPwUwHQAM1F45/gD1+tUdZ6qNqtqc2OjvbEbESVAUDgFH+cxAFX9nKo2qWqDqk5W1RvSP5B0ML+IPBc3wwaQRn6l3o0oIg0oBNXNqnoHAKjquh5//3MA96Q9DiJyE0jUe8/kDvOLyH8hZFiqky0REQA3AFimqtf0qDcVr4cAgE8B4F2mibISQNt0FphfRIEIIMPSPrN1JIAvAnhJRBYXa/8I4HMiMhOFJcPeAPDVJHZWJ462CIfhE+w1rq7OmG98905Tu675g6b28Ndnm9q/nmg7Od6YbdeBWrHRdm3MmGC7NnYfbn9Usy97yNTeWT5w7rvWJKvE0XuNMrX7V7ZGeu15h9vvyZ4jov1a/uCPq01t1bvbTe3C2ZMjbe/EaWNM7ZaFdh+uNe02b7froV33tO0E+vtZA19f7VoL7QO72TXYuhHt931g/q8rlpGq5hdQfoblKb8AfzKM+eXGDLPS7kZ8AnB+N+9Lc79EFI8kfAf5WsD8IgqH7xnGO8gT5Z34/66QiKikADKMky0iAsTvd4VERP3yPMM42SLKPfE+qIiISvM/wzjZIiKo50FFRNQf3zMsnMmWFD6Wjatx+GBTe/iaM03tY3/7n6bWtnWjqS37/R9N7eAWuybVC//+F6Y2dbRdf8pV63QsF+Vaf2rzGttpEqVLx7Wswbnf+rmp3fD9rwy4rVJca4G5umBeeHubqX3bsZ7XsMH2V/WbR9uOqb3H2O+ny+UPv2Zq/zxnH1O78cyDIm3vJ8+8ZWrPr7Jr0907fnivr0+ZYdeq22OUXVvu3pftunT7H2k7ocoi8P5dYfDKzC8gWoblKb8AfzIs7/kFMMOiCmeyRUQpqWAmQESUOf8zjJMtIoq0hAURkbc8zzBOtohyTkWgdYwCIgpTCBnm9+iIqDo8v96BiKhfnmcYJ1tEued/2zQRUWn+Z1g4ky3VXp0go4bYoXd0dZvaYU22e+KJFttB8vN/O9vUXN0trm6ZLW/b7pN/fHClqf3tR23nxT6O7hNX545Lx057HFEkva6Yi2uNty442pQczjzUrgXmWuNr3lOrTO17J+0baR/nHDbF1Bas2mJqx+49OtL2vjbLjvlvHd08v3mu9+/KfUvXmee4uqBm7G7XFRON9v2Mwve26eD1yS8g2QzLU34BfmdYnvILYIZFFc5ki4jS43lQERH1y/MMq2h0InJCUgMhooyIxH/UAOYXUY0IIL8qPbN1A4CE7kpGRJnx/F1hSphfRLXC8wwbcLIlIneX+isAE5IdDhFlwffrHcrF/CLKB98zLMqZrY8D+AKAvlczCoBZiY+IiKpMvL8hYAWYX0Q1z/8MizLZegrAdlX9fd+/EJFXkh9SNEdMHmVqdRG73VwdPy6u9bdcXTDt79m1xW69baGpPb3Mrg316pPPmdoJf3mMqbnWuKobZNdM6+5oN7UszJlqu2DuX9la9vbOarbdNydPHxvptT/4o+0EmjHRdsas2rTd1KJ280T93etr284OUxvjWAvvCwc3lbX9SAJYV6wCXuYXkH6GMb8qk2SG1Wp+AcywqAacbKnqyf383VHJDoeIqs//e9SUi/lFlAf+Z1hioxORJ5PaFhFVmdTFe9QY5hdR4DzPryTvszU0wW0RURX5fnFpFTC/iALme4YlOdlK7lawRFQ94v8p+CpgfhGFKoAMS3V0IjJFRBaIyDIRWSoiFxTr40XkIRFZXvzvuDTHQUQDSPimgCJykoi8IiIrROSSKhxB4phfRAFJ+KamSWdYkme2XEfQCeAiVV0kIqMAPCciDwE4B8Ajqnp18SAuAXBxgmN5X2u77dp5Ya1dQ+rrs223yOfOPc3UfnXDb03N1eHz3jtvmtqyx2zNZdFza2zR0c3z5v/8H1M7/LuPmdrqF18YcJ+utdWienWT7SDab5ztRnGZuYftqmkaUf6v5bfuW25q7Y7OraVv2Q6srx45LdI+Vrba7pt177WZmqsDqe/6eLcusmuN/fi0/SONIznJvisUkXoA1wM4AUALgGdF5G5V/VNiO0mel/kFRMuwvOcXkE2G5T2/AGZYVLFHJyKji+/sxovI+B5/9cW+z1XVtaq6qPjnrQCWAZgE4HQA84tPmw/gjNgjJ6LEqNTFegxgFoAVqvqaqrYDuAWFf/OZY34R1aYE8wtIIcMiT8FF5KsA/hnADvz5+gYFsA8AqOqSAV4/FcAhAJ4GsLuqri2+bq2I7FbiNXMBzAWAKVPsLJuIEhL/XWGjiPS8GdM8VZ1X/PMkAD1vDtQC4KMVjK5izC+iGhcvw/rLLyCFDItzvvNbAD6oqhvi7kRERgK4HcCFqrpFIn5mWjz4eQBw6KGH8gJWohSoCDT+4qwbVLW5xN+5Npb1v1/mF1GNKiPD+ssvIIUMizMVXAnA3qJ2ACLSgEJQ3ayqdxTL60Skqfj3TQDsrYmJqDoU0JiPAbQA6HkqZzIAx4U8VcX8IqpVyeYXkEKGxTmzdSmAP4rI0wDev5pOVb9e6gVSeAt4A4BlqnpNj7+6G8DZAK4u/veuOIMmoiQpuiMmUETPApghItMAvAXgLACfT3IHZWB+EdUs/zMszmTr/wF4FMBLAKItLggcicKFpy+JyOJi7R9RCKnfiMi5AN4EcGaMcQAAlm3YaWofbBxiaiMb7Mm7C2dPdmzR/qD2bhxuanM+e0K0AVZg8KDyuyqe/PYxjmrv2vJW23mzZovtRnH98u7stD/6Nzfbn4Wrk8e1FtiidfZkQyXdPN8/ZYapfeX2paZ25qH2d2DjDtulU8l6jlFcedIBZb+223mmuzxJxpSqdorI+QAeAFAP4EZVtT+E6vIqv4CkMyw/+QX4k2F5zy+AGRZVnN+KTlX9ZpyNq+oTcH/2CQDHx9kWEaVDATg6uivbpup9AO5LdqsVYX4R1agQMizOZGtBsbvmt+h9Gn5jUoMhomxosqfgfcT8IqphvmdYnMnWrs8rL+1Re791mojClMa7Qg8xv4hqVAgZFnmyparRbk9LRMHxPKcqxvwiqm2+Z9iAky0ROU5VHxWRT7v+vkc7NBGFSP1/V1gu5hdRDgSQYVHObB2FQhfPaShMHqXPfzMJK1fXjsugiM0OnY4f1NzmSaa28PX0L/E49gB7Q+qn17xnaq7uE5e+HTQzxtouG1ctqgWr7FqT63fY9dZ2G1Zvah/ezXZMuTz4ul0LzOXEaWNM7Wsft58Ubdhuv3cdjjXI0jZlVJLLk5ZHAXR5fr1DBbzMLyDZDMtTfgH+ZFje8wtghkUV5bu0VUS+CWAJ/hxSgP9n7YgoIs9zqhLML6Ic8D3Doky2di1rvj+Aw1C4gZ+g8E7x8ZTGRURV5Psp+Aowv4hywPcMG3CypapXAoCIPAjgUFXdWvz6CgC3pjo6IkpdYQkLz5OqTMwvotoXQobF+bB1LwA9b9vbDmBqoqMhokxkc7VHVTG/iGqY7xkWZ7L13wCeEZE7Ubje4VMA5qcyKiKqKs/fFCaB+UVUw3zPsDj32bpKRO4H8PFi6Uuq+nw6w7IEA3fluLpHGodFPUT7k9rcZrdXDY3DbVdN1M4dl9b23nP+Jevtel4udY7v9xGTRprasXuPLmtcQPRu0a6IH8i71gJzdTM9/Eb538+0udYLu/6ZFlP7+1m226wchRsCep5UFQohv4BKMiw/+QX4k2HMLzdmmBWrZ1NVFwFYlNJYiCgjfsdUMphfRLXL9wzL/gYZRJQ53zt5iIj643uGcbJFRN5f70BE1B/fM4yTLaKcUyi6vT8JT0TkFkKGcbJFlHfq/7tCIqKSAsiwYCdbrrXAxg+1a++5OhTqJFoLyb89ujL2uOI6q3mKqUXtXInqydV23a++Rg2xvwqHTx5laovetmuc7TXarvG2YUenrTnW8/rYZNsZ5DJ9vF2DbOXGaB1JGvHnnQVX91GH4+f/4mr7PCTUyQP4f71DLUo7w/KUX0A2GZb3/AKYYVEFO9kiomQo/H9XSERUSggZxskWEXl/vQMRUX98zzBOtojI+3eFRET98dANnwwAABeUSURBVD3DONkiyrkQ7r5MRFRKCBlWl+bGReRGEVkvIkt61K4QkbdEZHHxcUqaYyCiASjQ1R3vkRfMMKIAxMywLKR9ZusmANcB+GWf+rWq+v1KNuzqxnHNbKPOdZ9baztUzvhQk6mNGmK7ha57/DVTO/PQyaZW71qoyyHpbp4otnfYddR+t2JTpNdudHTtuH4WruNydbI01Nv3AMdPG2Nqb7TuiLSPqMfhs887ur6SEsK7wgzdhAAyLO/5BWSTYcyv6PKeYalOtlT1cRGZmuY+iKhSii7PgyorzDCiEPifYal+jNiP80XkxeIp+nGlniQic0VkoYgsfGfDhmqOjyg3dr0rjPOohIicKSJLRaRbRJqTOYqqGzDDmF9E1RE3wypRbn5lMdn6KYDpAGYCWAvgB6WeqKrzVLVZVZsnNjZWa3xE+VL9a7aWAPg0gMcr3lI2ImUY84uoSqp7zVZZ+VX1bkRVXbfrzyLycwD3VHsMRPRn1b7eQVWXAYB4fmfsUphhRH6pZoaVm19Vn2yJSJOqri1++SkUZolElCHfr3fwCTOMyD++Z1iqky0R+RWAYwA0ikgLgO8AOEZEZqIwGX0DwFejbEtFeq0R5ZxTOr7XW9rsOcMxQ+ynpx9pGhFlGBXdo9bVfTRn6mhTc3W4pK2SDqKOhHtpXdtr2WLXVWwaadcza9myM9mxVPB9aXB0b0XZXr3jl/vYve3vSVIK7wpjv6xRRBb2+Hqeqs7b9YWIPAxgD8frLlPVu2LvLSNJZVjf/AKSzbC85xfgT4blPb+AIDKs6vmVdjfi5xzlG9LcJxHFpGX9j2qDqpa8OFRV51Q2KD8ww4gCED/Dqp5fvIM8Uc4pKu/QISLKSggZltWtH4jII10a71EJEflU8SO5wwHcKyIPJHEMRJRfvucXz2wR5Zxq8tfd9b8/vRPAnVXbIRHVtGpmWLn5xckWUc6FsNQFEVEpIWRYMJMtUYWU8c10dR7+8yN2LbDLj98n0vbqHR05I4c2mNodi98ytUuO3y/SPk6ePtbUVrbabpbXHWtrVfMMRVpca7Ct2dpmapvb7HpmlXB12rg6cqK+NkrnTtSOn3sd66O51mArV6Wn1ql/5eYXEC3DmF/+yFN+ldoHM8wKZrJFROkI4V0hEVEpIWQYJ1tEeaeK7grux0NElKkAMoyTLaKcU/h/Cp6IqJQQMoyTLSLy/hQ8EVF/fM8wTraIcq7wrtDvoCIiKiWEDAtmsqUAOnt8LwfFW3C7l6idO+592B+oq/ukacwwU5s+1nb99F0vDYCza2mfcYMjbe/VTe2mtnLjdlPzmWvZhY07bDdTVlw/76jK7VB0qWQtuF4U3l/vELq++QWkn2HMr2zkPb8AZphLMJMtIkpHCNc7EBGVEkKGcbJFRN5f70BE1B/fM4yTLaKcU6j31zsQEZUSQoZxskWUd5rgtRNERNUWQIZxskWUcwr/g4qIqJQQMiyYydb2jm48//afu1I+0jTCPMfVAOH6HLfO0UFT5+jScVm83q7ndcmx001t4nDXt9buw3Xms4ImJezn6Pp5a8vOXl+7Okp8/0X1yZgh9mfrWufM9XvWUN/760rWgqukq6gnDeBdYej65heQbIblKb8AZlglkswvgBkWVTCTLSJKj+9BRUTUH98zjJMtopxTqPdBRURUSggZxskWUd4FcAqeiKikADKMky2inAvh4lIiolJCyLC6NDcuIjeKyHoRWdKjNl5EHhKR5cX/jktzDETUv10Xl8Z55AUzjMh/cTMsC2mf2boJwHUAftmjdgmAR1T1ahG5pPj1xQNtaHhDHQ7ZY/j7X1dyt1h3d4/jeY6+mqljh5ja6MG2RcPRyOFq5qlI1HXJjpk6ZsDnLFi1xdR2dpbfZeK7qGt8uZ7n6tw5cdoYU1u6oc3Upo213Vbmde/YjrGObvuzGD/Mri1XrjxNoGK6CQlkWN/8ApLNsDzlV6nn5SnDfM4vgBnmkuqZLVV9HMDGPuXTAcwv/nk+gDPSHAMR9W/XxaW+vzPMAjOMyH9xMywLWVyztbuqrgUAVV0rIrtlMAYiKlIFOnM0gUoAM4zIIyFkmNcXyIvIXABzAWDKlCkZj4aoduXpbFW1ML+Iqsf3DEv1Y8QS1olIEwAU/7u+1BNVdZ6qNqtqc2NjY9UGSJQnvEA+tkgZxvwiqo4QLpDPYrJ1N4Czi38+G8BdGYyBiHroUo31yDlmGJFnfM+vVD9GFJFfATgGQKOItAD4DoCrAfxGRM4F8CaAMyNurddaTVHXPPyXR18ztZZN203t53/5QVNz7aO9y9Zcz6t3tvNYK1ttx8cMR8fHwjXvmZqrk2P6WFtzde70FbW7hYBbFq42Ndd6Yx/d0659F4VrzbyIv05lqfbdl0XkPwCcBqAdwEoAX1LV1qoNIIbkMkxMPiWZYXnPL4AZFlXa+QXUdoaVm1+pTrZU9XMl/ur4NPdLRNFlsIjrQwAuVdVOEfkegEsR4fYvWWCGEfmvyhlWVn55fYE8EaVPAbR3dlVvf6oP9vjyKQCfqdrOiajmVDPDys0vTraI8k4zvej9ywB+ndXOiagGZJdhkfOLky2inCtzXbFGEVnY4+t5qjpv1xci8jCAPRyvu0xV7yo+5zIAnQBujrtzIqJdysiwqucXJ1tEOVfmDQE3qGpz6W3qnP5eLCJnAzgVwPGqbG8kovKVkWFVz69gJlsCRV2Pxbnca4FZh08bb2q3OroRXQY59jF+qF1H7J7lfVfzAE6dYffr2p6rc8dlmmNNs/HDXD++gX/uT7RsM7VxQ20X0GFNw00tqvtXJttcdvL0sab24OubTa0ap5Jd+3imxR5vud08Ubt2onZpDaTMM1tlE5GTULig9GhVjfaPMXB98wtIP8NqNb+A8DIsT/kF1HaGlZtfwUy2iCgl1e9GvA7AEAAPSSGVn1LV86o5ACKqIdXNsLLyi5Mtopyr9n22VHXfqu2MiGpeNTOs3PziZIuIuAQPEQXN9wzjZIso5zK4qSkRUWJCyDBOtogI6nlQERH1x/cMq6nJlqsDYs402wWy4l27TpdLN+wGb1/2jqk17znG1FwdPp/cb4Kp9e1QKmWso4soqnd39u5x2m2E7SDab1y0riKXxet3lP3aSgxvsN+TrW2dkV47dJBdg32vMUMjvXZ7u71T8ebtHaa2vLXd1Bat2dLr67/6QKN5TlIdOlGpAt2eB1VeJJlhtZpfQG1kWC3kF8AMi6qmJltEVA4Fb3VFROHyP8M42SIi70/BExH1x/cM42SLKO8COAVPRFRSABnGyRZRzikAdd26nIgoACFkGCdbROT99Q5ERP3xPcOCnWy5Om1cyzG5umX+rnnPaPtw/PCOmGI7g5pG2m/j3mPs2mKu7bkm4641yOqiLjbl0PeXcN9xdp0y15pk6tinq8tkg6OTJWnrttsOmiOnjDK1SrpgXt7YVvZrJ46239NfPPWmqW18r/c+Wnfa791+E+yaZK7OpYb68n8negngFHwtSjvDajW/gPAyrFbzC2CGRRXsZIuIkqLeX1xKRFSa/xnGyRZRzhWud/A7qIiISgkhwzjZIso7dX9EREQUhAAyjJMtIvL+XSERUX98z7DMJlsi8gaArQC6AHSqanNWYyHKO9+DykfMMCJ/+J5hWZ/ZOlZVN5TzQleHzvwX1pta0yjbZXHiNLsWWFQvrttq9zFynKm5um9cpzmj/npEXYPMpXFY3y6QaNtydcZ8677lpnb8/hPLGVYsi9ba73t9nf0eR/3Zvra5/O4j134fXrqurG01jbLrmW11rF3mqiVFVb3v5PFYUBlWG/kVfc++ZFit5hfADIsq68kWEXnA93vUEBH1x/cMs0uHV48CeFBEnhORuRmOgyj3tDvegwAww4i84Xt+ZXlm60hVXSMiuwF4SEReVtXHez6hGGBzAWDKlClZjJGo5mkANwT0VL8Zxvwiqo4QMiyzM1uquqb43/UA7gQwy/GceararKrNExsbqz1EotzQbo31oIEzjPlFVD2+51cmky0RGSEio3b9GcCJAJZkMRai3FNOtuJihhF5JGaGZSGrjxF3B3CnFDpeBgH4H1X9XZwNuNYVW7dlZ6TanGl2fbCo3TL3Obo2Tt7XriMW1b3LN5rakjVbTO3Lh002tT1HlPfjc33vXJ1Grq6Vq0+eYWoPvdYaab8N9XZuP3SQre02YrCprd1q1+Ta3mG7W5ZusM/7YKPt5mpwHJuqHcsgRyPU8MG2uHVnp31iypJbVkzR1cULsWLyKsPylF+APxnG/KpMnjIsk8mWqr4G4OAs9k1Efaj/96jxDTOMyCMBZBhv/UBE3l9cSkTUH98zjJMtIvL+HjVERP3xPcM42SLKOdXqXjQqIv8C4HQA3QDWAzhnV2cfEVFc1cywcvMry5uaEpEnurs11qNC/6GqH1bVmQDuAXB55UdARHnme34Fc2ZLIb06UFa22q6N8SNtF8iHdhtlaq6uFUdzBwY5ateffoCpdVZw+vIjTaNN7YkV75raU6s3m9qnD5hgavNftGur9e2OOWEfRzemYy20li12/a0f/e/rpuZaV+zk6XYfLq9uaje1be22M+bovezP0dW5M9j1Q3OYNKrB1FwdXctb7fjaOv3oejlw4sjEtqXd6a1bZval2rNdbQSiL7EXrL75BSSbYXnKL8CfDGN+VSbEDCs3v4KZbBFRSlSrOtkCABG5CsDfANgM4Niq7pyIakuVM6yc/OLHiEQ5pygEVZwHgEYRWdjj0WttQBF5WESWOB6nA4CqXqaqUwDcDOD86h81EdWKuBmGDPKLZ7aI8k4B7Yr9rnCDqjaX3KTqnIjb+R8A9wL4TtwBEBEBKCfDqp5fnGwR5V7VT8HPUNXlxS8/CeDlqu2ciGpQ9TKs3PziZIso76p/zdbVIrI/Cq3TqwCcV82dE1GNqW6GlZVfwU626h2dJ2OG2A6Nw/a03Q6Ol+InT7eYWpejRfTC2XaNL1cDyRMt20xt9iQ7lqaR9kdw0dHTTG3MEMciVw7Txw+PtN++XB2aY4fafbrWFbvsgRWmdtTeY0xthOMbNWO8XffrD6ttB43LsAZ7yeEaxxpkM8babqaoa8m5tHvSzfPudttpVa4qdyP+ZdV25rEkMyzv+QVkk2HMr8qEmGHl5lewky0iSoZW+WNEIqIkhZBhnGwR5Z1W98wWEVGiAsgwTraIck/R7XlQERGV5n+GcbJFlHcZ3NSUiCgxAWQYJ1tEOafw/xQ8EVEpIWRYOJMt6d2BM3WM7dCYOqYx0qZcS4F9bZbt0nFxrSPm6ub52GTbQdMZsYGkaYT9sUR97ZB62+HS96XRVt8ChkVcp+vI6XaNs3Xv2fXB9hljO63E8f3s6LLdMuro3HJt7+1ttpvHpc3RkLOt3Ra/+8ArpjZ8sO1w2t5e3j/0/35qlal9cfbekV67uc1+j8uiWs5NTSkOsR2ESWZYnvIL8CfD8p5fADMsqnAmW0SUGt/fFRIR9cf3DONkiyjvArjegYiopAAyjJMtotzzP6iIiErzP8M42SLKucLFpf7cVZqIKI4QMsxejVglInKSiLwiIitE5JKsxkGUe8VT8HEexAwj8kbMDMtCJme2RKQewPUATgDQAuBZEblbVf8UdRuutaFcHR+uzkP3mKI9z7WemWMJMuf4XI0xri6dDkctavfNIXvYtcWSVOc4/lP2HWdqz7+93dSmjLbdNw2OAzt279Gm5voeu74nR0RcR+3sX71gavV1douDHLW07TfB/gxd3/fhjrXVysUJVDy+ZRjzK7okM4z55cYMs7L6GHEWgBWq+hoAiMgtAE4HEDmoiCgh6v/dlz3EDCPyRQAZltVkaxKA1T2+bgHw0YzGQpRrCnh/jxoPMcOIPBFChmU12XKd2zQnWkVkLoC5ADBlypS0x0SUT9qN7s72rEcRmgEzjPlFVCUBZFhWF8i3AOiZPpMBrOn7JFWdp6rNqtrcOHFi1QZHlDe8QD62ATOM+UVUPb7nV1Zntp4FMENEpgF4C8BZAD6f0ViI8i2AGwJ6iBlG5IsAMkw0arte0jsWOQXADwHUA7hRVa8a4PnvAFgFoBHAhvRHmLpaOI5aOAagNo6j5zHsraqRT6WIyO+Kr49jg6qeFPM1NSVOhvXIL6D2ft9CVgvHUQvHAPz5OGLlF1BWhlU9vzKbbJVLRBaqanPW46hULRxHLRwDUBvHUQvHkBe18LOqhWMAauM4auEYgNo5jlIyu6kpERERUR5wskVERESUohAnW/OyHkBCauE4auEYgNo4jlo4hryohZ9VLRwDUBvHUQvHANTOcTgFd80WERERUUhCPLNFREREFIygJlsicpKIvCIiK0TkkqzHE5WI3Cgi60VkSY/aeBF5SESWF/9rV0L1iIhMEZEFIrJMRJaKyAXFejDHISJDReQZEXmheAxXFuvTROTp4jH8WkQGZz3WgYhIvYg8LyL3FL8O7hjyhvmVnVrIL4AZFrJgJlsiUg/gegAnA/gAgM+JyAeyHVVkNwHoe0+PSwA8oqozADxS/NpnnQAuUtUDAcwG8LXi9z+k42gDcJyqHgxgJoCTRGQ2gO8BuLZ4DJsAnJvhGKO6AMCyHl+HeAy5wfzKXC3kF8AMC1Ywky0AswCsUNXXVLUdwC0ATs94TJGo6uMANvYpnw5gfvHP8wGcUdVBxaSqa1V1UfHPW1H4RzIJAR2HFmwrftlQfCiA4wDcVqx7fQwAICKTAXwCwC+KXwsCO4YcYn5lqBbyC2CGhSykydYkAKt7fN1SrIVqd1VdCxSCAMBuGY8nMhGZCuAQAE8jsOMonrpeDGA9gIcArATQqqqdxaeE8Hv1QwD/AKC7+PUEhHcMecP88kTI+QUww0IV0mRLHDW2UlaZiIwEcDuAC1V1S9bjiUtVu1R1JgoLB88CcKDradUdVXQiciqA9ar6XM+y46neHkNO8WfkgdDzC2CGhSqrhajL0QJgSo+vJwNYk9FYkrBORJpUda2INKHwLsVrItKAQlDdrKp3FMvBHQcAqGqriDyGwvUbY0VkUPFdle+/V0cC+GRxXb6hAEaj8C4xpGPII+ZXxmopvwBmWGhCOrP1LIAZxY6FwQDOAnB3xmOqxN0Azi7++WwAd2U4lgEVP1O/AcAyVb2mx18FcxwiMlFExhb/PAzAHBSu3VgA4DPFp3l9DKp6qapOVtWpKPwbeFRV/xoBHUNOMb8yVAv5BTDDgqaqwTwAnALgVRQ+o74s6/HEGPevAKwF0IHCO9xzUfiM+hEAy4v/HZ/1OAc4ho+hcFr3RQCLi49TQjoOAB8G8HzxGJYAuLxY3wfAMwBWALgVwJCsxxrxeI4BcE/Ix5CnB/Mr02MIPr+Kx8EMC/TBO8gTERERpSikjxGJiIiIgsPJFhEREVGKONkiIiIiShEnW0REREQp4mSLiIiIKEWcbBERERGlKKQ7yFOViMgVALahcGffx1X14Qq2NRPAT4vb6gJwlar+OolxEhH1xfwiH3GyRSWp6uUJbGY7gL9R1eUisieA50TkAVVtTWDbREROzC/yCT9GJACAiFwmIq+IyMMA9i/WbhKRzxT//IaI/KuIPCkiC0XkUBF5QERWish5pbarqq+q6vLin9egsPbYxCocEhHlBPOLfMczWwQR+QgKa1QdgsLvxCIAzzmeulpVDxeRawHchMKCokMBLAXwswj7mQVgMArLlRARVYz5RSHgZIsA4OMA7lTV7QAgIqUWyN1VfwnASFXdCmCriOwUkbH9nVoXkSYA/w3gbFXtTnDsRJRvzC/yHj9GpF2iLJLZVvxvd48/7/q65MRdREYDuBfAt1X1qbJHSETkxvwir3GyRQDwOIBPicgwERkF4LSkNiwigwHcCeCXqnprUtslIipifpH3+DEiQVUXicivASwGsArA/ya4+c8COArABBE5p1g7R1UXJ7gPIsop5heFQFSjnH0lIiIionLwY0QiIiKiFPFjREqEiHwIhW6dntpU9aNZjIeIKCrmF6WNHyMSERERpYgfIxIRERGliJMtIiIiohRxskVERESUIk62iIiIiFLEyRYRERFRiv4/Ykg0g7t8HJcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=[10, 10])\n",
    "ax1 = plt.subplot(221)\n",
    "ax2 = plt.subplot(222)\n",
    "xds1[-1,:,:].plot(ax=ax1)\n",
    "xds2[-1,:,:].plot(ax=ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read/Write Validation\n",
    "\n",
    "Let's make sure that we can read a .pfb file, write the data we read back to disk, then read that data back and make sure it's the same! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "pfb_c_outfile = 'test_pfio_write.pfb'\n",
    "pfb_f_outfile = 'test_fort_io_write.pfb'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "pfio.pfwrite(np.flip(pfb_c_data, axis=1), pfb_c_outfile, 0,0,0,1000,1000,1000)\n",
    "pf_fort_io.pfb_write(np.transpose(pfb_f_data, (2,1,0)), pfb_f_outfile, 0,0,0,1000,1000,1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "pfb_c_data_read = pfio.pfread(pfb_c_outfile)\n",
    "pfb_f_data_read = np.zeros((41,41,50), order='F')\n",
    "pf_fort_io.pfb_read(pfb_f_data_read, pfb_f_outfile)\n",
    "\n",
    "pfb_c_data_read = np.flip(pfb_c_data_read, axis=1)\n",
    "pfb_f_data_read = np.transpose(pfb_f_data_read, (2,1,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.remove(pfb_c_outfile)\n",
    "os.remove(pfb_f_outfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "xds_c_read = xr.DataArray(pfb_c_data_read)\n",
    "xds_f_read = xr.DataArray(pfb_f_data_read)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f2eadf50250>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAEmCAYAAABRZIXOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZQc5Xkv/u/TPT2bZkajmdEyaEECJEBswsgyhIAxeAFiszh2rnHi4IQb2b8TJ3ac3BgnObGdxDn2vXbwPTe/OJEvNtg3Pww29oGLgzEBDDExwkIItCGEFrSNds1Io9l6up/fH90SM/O8Ule3qqbf6vp+zumD+pmarqruni9vV9dTr6gqiIiIiChcqWpvABEREVEt4iCLiIiIKAIcZBERERFFgIMsIiIioghwkEVEREQUAQ6yiIiIiCIwKYMsEUmLyMsi8ljx/gIRWSkim0XkQRGpn4ztICIqF/OLiColk3GdLBH5LIClANpU9f0i8hCAH6nq90XknwG8oqrfPN1j1E9t0sZZbSfvL2rsq3h7XHu8L9dgau2pQVNrlHDHpaPIm9qebKupNadGTK0rna1onTnHOl873mVq0pc2tbx9mrCoY6+p1Yv9XdfzfjRvn88Rtb/bIKOm1pqy+yEQx3rDfo/bdfTm7TanHc9zayrItthlsmofawR2nQCwae3wQVWdHmBFAID3vWuKHjqcC7o4AOClV4efUNUby/qlmIoiv4BwMyxJ+QX4k2HMr1OpPMPKzS+g/AybzPyqi3oFIjIHwG8A+DKAz4qIALgewEeLi9wP4IsAThtSjbPa8I5/+ejJ+09c8JOKtynvePP8j8OLTO3W1ldM7YLMlIrX63Iof9zUvrD3Xab2tpY3Te332/ZVtM6jeRu+v/biXabW+NOp9nfPtX88j97xP0xtXp0NWtfz/vhAs6ntyU4ztXPq95vatY3DppYR+5bOqg24M5F2/I/q4X67ze1p+9re0GT/ZzNRzhFGPTn7mu0cdb8Xr12wxb5ZTuPg4RxWPjGnnF9BpnuL/T9aDYoqv4BwMyxJ+QX4k2HML7czybBy8wsoP8MmM78iH2QB+AaAPwdw4h3bCaBX9eQ7ZxeA2ZOwHUTkpM5QJADML6IY8DfDIj0nS0TeD2C/qr40tuxY1HkMUkSWi8gqEVmV7bOjYCI6cwogDy3rlgTML6J4KDfDJlPUR7KuBnCLiNwMoBFAGwqfDNtFpK74aXAOgD2uX1bVFQBWAEDb+TOTkexEVeD6KpeYX0Rx4WuGRXokS1U/r6pzVHU+gI8AeFpVfxvAMwA+VFzsTgCPRLkdRHRqCkVOy7slAfOLKB7KzbDJNBnnZLl8DsD3ReTvALwM4N5SvyAC1I3pxvinXnuS27OHzze1EUfXxMBoxtRu715jaosyTabmGi2/e/1tprbtzRmmlm6yJzBmGmzN9X3E4sX2w7JrW1IBxs3PDnWaWs7RJTNsz4XEaLN9g97Xu8zU5tUfMrW+nD1J9DtvXGlqN87baGofaNlsammxj+d6Tl7L2i6m3rx9ba9uDPZJyPXd/z/vuNbU3j1zk6nd0PRaycd3nZg6p851wrI9MbVSSfkKMCRnnF9AuBmWpPwC/Mkw5pcbM+wtkzbIUtWfA/h58d9bAdh3NRFNOgWQ8zSgfMH8IvKXzxlWrSNZROQRXz8FEhEF4WuGcZBFlHAKJOY8KyKqPT5nGAdZRORpXw4RUTC+ZhgHWUQJp1Bvz2cgIirF5wyLzSArIznMaOw/ef9/bbjOLDPcY7sXGg/YLoehc+x0BnPn2U4Sl28cOdfUdh1qN7XmbbaDMZW1NVcrTt3Vh03tmmbbnQLYSbhcUzCMYvycThdk+s0ybU1Dpnao29EN0mznh3q5d66pXX6WnRlhSO3+33z2BlN7V6utzUq3mJqrE+cnA3a51pTtzrq6wTVvmn2v7Bi13S//88B1prZtt51qa+dUO5WGXaNdZ9Drvbi7dSqgQM7PfKoZE/MLCDfDkpRfgD8ZlvT8KqyVGXY6sRlkEVE0CldLJiKKJ58zjIMsosQT5JxXNyIiigN/M4yDLKKEUwB5Tw+1ExGV4nOGcZBFRN5+CiQiCsLXDOMgiyjhCldL9jOgiIhK8TnDYjPIakqPYEnrzpP3V+Jss8yUXbbLIWUbCVE/baDi7VjcuNvUzp5uu2n2iO2aaOizj3f4HSOmdkWH7RSamrJdNynY+atc77M/73nHuPuZlD1F8Ohgo6k1nX3M1P72YjsX7i1T7HLujhPbEQTsdNQs10mNG7O2m2jT0EJTe1eLnUssLcHe+jtHW01t/7CtpTN2C4+M2PfAsI7vCmoS22EVVNB53oLIq58BVSsm5hdQnQyrhfwCfMqwZOcXwAwrJTaDLCKKhs+fAomISvE5w8IbRhJRLCkEOaTKupUiIp8WkXUisl5EPlOsfVFEdovImuLt5sh3johqXrkZFkRYGcYjWUQU6qF2EbkYwB8AWAZgBMBPReQnxR/fo6pfC21lRETwN8M4yCJKuAgOtV8I4AVVHQAAEXkWwO1hroCI6ASfM4xfFxIlnEKQ1bqybiWsA3CtiHSKSDOAmwGcmLfkUyLyqoh8W0SmRblfRJQM5WZYAKFlWGyOZB3PNeCF3nNO3m9rtp0ZhzvsvE+prB3dXjXTdtjc1Gy7dVydD+91zI/14vQtpvbQqJ0LK+uYpmnGLNuyc0On7SZxzfEUdH6o+Y3ju32+ueEas8zIYduZM+e8HlN7X7PtREo55iALum27HHNrNYp9zTZnm01tz2i3qW0f6jS19jbbAZVCvam5tvnyevt6XzPNzsM2NTNoauc37zO1rI5fR9MZfPgK+hwHUcGnwC4RWTXm/gpVXQEAqrpRRL4K4EkA/QBeATAK4JsA/haFD55/C+DrAH7/DDc9FibmFxBuhiUpvwB/Mizp+QXENsNOmV9AuBkWm0EWEUVDVZDTsg9qH1TVpad+TL0XwL0AICJ/D2CXqp5MbRH5FoDHKthcIqJxKsiw0+ZX4THDyTAOsogI+ZDbn0VkhqruF5F5AD4I4CoR6VbVE4cWbkfhkDwR0RnzNcM4yCJKuMJJo6GfnvmwiHQCyAL4Q1U9IiLfE5ElxVVuB/CJsFdKRMnjc4ZxkEWUeBV9XXhaqmpOmlHVj4W6EiIiAD5nWKSDLBFpBPAcgIbiun6oql8QkfsAvBPAibMmP66qa6LcFiJyUwB5NhobzC+iePA5w6I+kjUM4HpV7ReRDIBfiMjjxZ/9N1X9YdAHGjjWiJd+sejk/VyTmmWaeu13skMX266J/zrj2UDr3D5q57Pakm03NdcIerjLbl9q2G7fpdMOmNrNU15zbI3tznF1tsyrs/NSXdw4fn6tfM5uh6bs9jbX2XmqNmXtcs8en29qvzhynqld0LLX1B5+Y4mpNdTb9V4/x9ERU2df2xn19jVbUGc7e1y2ZG2H1sx02tQW1tv9uKjTdqwurrevz8R5vlzdNTlHB8+hvO0Syji6mCqV83TeryqLLL+A6DOsVvML8CfDkp5fADOslEgHWaqqwMkZNTPFm32HE1HVnJiSgsZjfhHFg88ZFvlWiUhaRNYA2A/gSVVdWfzRl4sX9LpH5Aym8SaiM5bXVFm3pGB+EcWDr/kV+dpUNaeqSwDMAbCsOCfQ5wFcAODtADoAfM71uyKyXERWiciq3HF76JKIztyJzpywJ1itBcwvIv+Vm2GTadLWpqq9AH4O4EZV7dGCYQDfQWESRtfvrFDVpaq6ND3FcblhIjpjCkFOy7slDfOLyF/lZthkinSQJSLTRaS9+O8mAO8G8JqIdBdrAuA28KKERFWVR6qsWxIwv4jiw9f8irq7sBvA/SKSRmFA95CqPiYiT4vIdAACYA2AT5Z6oLrjwKwX3upYyDXY0Wiu3p6T2nat7da4pH7Y1F4atk/FP+272dR6R5pMbf+AnW8sPWS3b+QsOwfV78543tRc83y5TE/b+atc3tU0vtvl60ttU9T399sP45sOzzC1v9xuJyIfzGZMbc+Rqab2WsY+3uARO9/Y8IDd/0eOXWpqdfU5U/vgwldM7YFjXaY2o+6oqeXQYWrPD9purLkZO5fa1Y2uObjse6U/P/69t2nEvu/WD882tX2jbaZ2KGs7sQoeOkXdTdXdYUbR5RcQboYlKb8AfzIs6fkFRJFh5eUX4HeGRd1d+CqAyx3166NcLxGVQ0KfkqIWML+I4sLfDOMV34kSTuHvp0AiolJ8zjAOsogoUR2DRFR7fM0wDrKIEk4hyCewY5CIaoPPGcZBFhF5+ymQiCgIXzMsNoMsrQOGpr31JIptzMBImx3JtqRcXRPWS0PzTW1Vz1xTGzhmL+4sB22tzjH5RnO7natqVtp2DhXmoy3tcN52+2wYsZ0t52R6x91/eeAis8wre88ytaEB2/0zrcnOj3XHnF+ZWnqufd5XHZtvauct2m9qrx+fZWo/W7fY1EZNBbioaZepPX90oanNarDdOX/Vtd7UHle73GxHZ09f3r73pqZs51FLavxr25qy74nOun5Ty4jd28uadpgaAHzdWT01BRJ1FfdqmJhfQLgZlqT8AvzJsKTnFxB+hpWbX4DfGRabQRYRRUWQ87Qzh4ioNH8zjIMsooTz+VMgEVEpPmcYB1lE5O2nQCKiIHzNMA6yiBJOVbz9FEhEVIrPGcZBFhF5eyE/IqIgfM0wDrKIEk4Bb6ekICIqxecMi80gS9tzGL39yMn7A4O2NbexwU4kWp+2baO3v/ZbjuVsP/XxPjtBZscv7XpTWdvv3G/nx0RHi20dbk25Gnltm2zKcQ2QvTm73IahOaY2PT2+nXbjsW6zTH6NnQw1P8M+J3WzbVtzX67Z1P6sY4upfbytx9Scptnf/caUvab2+F7bxn1Bg13ujtkHTe2+o/Y52DV63NRusruGI46O+g0jtt15WYNdMCPj/+TOrbOv4aKMbX/OOVqx0xLWJzfx9lNgrZiYX0C4GZak/AI8z7AE5RfADCslNoMsIoqGAshqutqbQURUEZ8zjIMsooTzeUoKIqJSfM4wDrKICHlPp6QgIgrC1wzjIIso4VSBnKefAomISvE5wzjIIiJvD7UTEQXha4bFZpC1oOkgvnfpd07ezzo6CVydLnfvuM3U1j9pJ9x0zF2JujbbdZO1DTuok2Av7gXtdjLRmWnb7eMyqMOmtjPbYWov9J5jartHpo2731xnH2togZ2sNXUkY2qv75xpah0Ntuvoe5k+U/t57/mmdl37JlP7WKvtpvmMo2Pno21rTW39iO0wAmzHlqtL6PmhVlObU2c7bKal7Jvgsnr7/KWl9GsbtLvG1dUTlsL5DH4eaq8VE/MLCDfDkpRfgD8ZlvT8KizHDDud2AyyiCg6vk5JQUQUhK8ZxkEWUcIVJlf1M6CIiErxOcM4yCJKPH8PtRMRleZvhkW6VSLSKCIvisgrIrJeRL5UrC8QkZUisllEHhQJ+OUvEUUiDynrVoqIfFpE1hX/7j9TrHWIyJPFv/snRcSebOMR5hdRfISZX0B4GRb10G8YwPWqehmAJQBuFJErAXwVwD2quhDAEQB3RbwdRHQKJ9qfy7mdjohcDOAPACwDcBmA94vIQgB3A3iq+Hf/VPG+z5hfRDFQboaVEmaGRfp1oaoqgBOTGGWKNwVwPYCPFuv3A/gigG+e7rEaJYULM3Z+pbF6crYLY2jU7mK2xXbdpLL2iZf5tuMkvdh2tRw50GI3JmvHr0//x6WmdlHHhfZ37eYh3Wjn4EpnbC2/xW7LylnjO0c6uo7ZFQzaKQmm7LLPyXGx81Ttn23X+e+HF5va7KZeU3N1HKZgO2LucHTszErb9XY22tfnmUG7zTc02efu6kbHpF4BNaccc8KF+Bkm73hOwhTyofYLAbygqgMAICLPArgdwK0Arisucz+AnwP4XJgrDtNk5xdQeYYlKb8AfzKM+RVcUjMs8i8xRSQtImsA7AfwJIAtAHpV9UTT8S4AjulIiWgynJiSopwbgC4RWTXmtnzMQ64DcK2IdIpIM4CbAcwFMFNVewCg+N8Zk72v5WJ+Efmv3AzD6fMLCDHDIj/xXVVzAJaISDuAH6MwQjSLuX63uOPLAWDebJ6jTxSVoOcpjHFQVZe6fqCqG0XkqygMSvoBvALAcSUn/zG/iOKhzAw7ZX4B4WbYpJ2Or6q9KBxauxJAu8jJK5PNAbDnFL+zQlWXqurSrk4/Z9gmirsT7c9lfhI8/WOq3quqb1PVawEcBrAZwD4R6QaA4n/t1S09xfwi8le5GRboMUPKsKi7C6cXPwFCRJoAvBvARgDPAPhQcbE7ATwS5XYQ0enlNVXWrRQRmVH87zwAHwTwAIBHUfh7B2Lwd8/8IoqPMPMLCC/Doj6G3Q3gfhFJozCge0hVHxORDQC+LyJ/B+BlAPdGvB1EdCplfLorw8Mi0onCnCB/qKpHROQrAB4SkbsA7ADw4bBXGjLmF1EceJxhUXcXvgrgckd9KwqtkaH6z8E5ptY7bOdpynXYr1YXnbvb1D4792emtmHYnuO6cVa3qT3zpp0fMb3dzi1Vv9V2dTTvs10YmUE7+h6aZruVUqOuDo7xL/Nog50zbJZtagLUPtax823tmul2Xq4ZmaOmdmWTXW562q54xeErTe1/jbaZ2qenbTc11/xYTxy1XVGz6p43tYsyzaYWtokdNn35IbPMpqx9T6wdmmtq7WnbPVYJRUXnZJ3+MVWvcdQOAbgh1BVFaLLzC6g8w5KUX4A/GZb0/AKYYaXwbEwi8nZKCiKiIHzNMA6yiBLO53m/iIhK8TnDOMgiIm8DiogoCF8zjIMsooQ7cSE/IqI48jnDOMgiotBPGiUimky+ZlhNDbLWOTpzeg5MNbX6Vjs/1G/OWm1qF9TbufZeGZpnanMbj5jaeV12rqotV5gSBm0Jx1+3XTzNPcGu7ZGdYt9oZz0/fi2ZNw+YZQYW2w6jkan2AooNe+1bJue47sh1za+b2qKM7ZICbBfKF6avNbXXs/aZyqqd0ystdluub1tvakHmkStHpfNyTU3Z7WhP2X2dW3/I1I7lXM9nBdTfQ+1JU2mGJSm/AJ8yLNn5BTDDSqmpQRYRlc/nk0aJiErxOcM4yCIibwOKiCgIXzOMgyyihPP5pFEiolJ8zjAOsogI6mlAEREF4WuGcZBFlHCqwGjASVOJiHzjc4bFZpB1XBUvDr/VAdGbt/M09QzbLhxHswYymZypzazrM7VGxy+/c8omU9s6Mt3Ups6wczLtm2a3b93Rs0xtzWDG1Ppm2i4ZV5dktsc+L/m68fsxOrvTLNO7yK6z9xI7L1f7TNuxtOGY7erZ3WLn6lqUsY+XQrA/jAsyU0wtaEfMe5vs3FpB1+tax0B+xNRW9F1oaosb7XxyNzbZ353o3DrbdXQoZ/dhVr2dW61Svn4KrBUT8wsIN8OSlF+APxnG/HJjhr0lNoMsIoqKv+czEBGV5m+GcZBFRN5+CiQiCsLXDOMgiyjhfL7GDBFRKT5nGAdZREmnhRNHiYhiyeMM4yCLiLyd94uIKAhfMyw2g6xRTWPv6FvdLVm1m75v0M6ZlTts55bKzrDdhb0529XSJv2mdqltYEGj9JjaarULnl1v5wPLqu26mXmB7bjoGbSdPSN5+7uvDc8ytYOXjp9baspeu22jthkEbdOPm9ofLfq5qa3uP9vUto7MMLW07DW1qxpst0q/2trOUbuvs9KjpvajYxeY2vtbNpravDr7XsmqfTzXXGLNKfueOpK13UPtzbZDa+KfnGud67O21payH9Muqne8aBVQ+Hs+Q62YmF9AuBmWpPwC/MmwpOfXqdbLDHtLbAZZRBQVfztziIhK8zfDOMgiIm/PZyAiCsLXDOMgi4i8PdRORBSErxnGQRZRwqn6G1BERKX4nGGRTvYjInNF5BkR2Sgi60Xk08X6F0Vkt4isKd5ujnI7iOj08ipl3ZKA+UUUH77mV9RHskYB/KmqrhaRVgAviciTxZ/do6pfC/pAU1N53DLl2Mn7D/dPM8s019l5qjrm23mqjvbbbpX1g3NMbX/zm6bWnbZdGBdm7OPVY7tdr9qujve2rTW1RrH7sXvU7u+uETt/12jejpt3NI3/3f4dLWaZlGNKqtFDdl+P5203yO91/cLUvnvwalN7bdDODzZl2kpT2z1qu3p2Zu2+vti3wP5uv+1iun6Rna/NZX9+0NT+Y3Ceqbk6qq6Yss3UljWU/gzj6v55ov9iU3tfyzpTCzp/WRC+ns9QZZHlFxBuhiUpvwB/Mizp+QUww0qJdJClqj0Aeor/PiYiGwHMjnKdRFQ+Xw+1VxPziyg+fM2wSL8uHEtE5gO4HMCJof+nRORVEfm2iNiPKoXfWS4iq0Rk1cFD9tpWRHTmFALV8m5Jw/wi8le5GTaZJmWQJSItAB4G8BlVPQrgmwDOBbAEhU+KX3f9nqquUNWlqrq0q9Me4iSicGiZtyRhfhH5z9f8iry7UEQyKATUv6rqjwBAVfeN+fm3ADwW9XYQ0Sl43JlTbcwvohjwOMOi7i4UAPcC2Kiq/zCmPvbswdsB2DPiiGjy8FCWwfwiihFP8yvqI1lXA/gYgLUisqZY+wsAd4jIEhR2dzuAT5T7wIdytrvkotY9pvbODtuZMb3Ozq315kiXqT09YLs/FjfYdVxeb5/GhRm7fUcd3R9DeszUGsXR/VF/xNSeT9u5yV5oOMfUZswev9ye9jazzNZtM01NBux2pJE3Ndf+/0vezi02mLXzq12YseeqXNFgnxPA1trTdm6txk7b2bQo02Rqecd+uDqv1hy33Tm7h9pN7fK2Haa2sf6AqV2UGf8cuLpr/kvby6Y2NWVfC9c+VCrsT4Ei8icA/isKf+NrAfwegH8G8E4AfcXFPq6qa9yP4IXI8gsIN8OSlF+APxmW9PwCmGGlMizq7sJfAM6psf8tyvUSUXnCbH8WkdkA/hjAYlUdFJGHAHyk+OP/pqo/DG9t0WF+EcWHrxnGK74TJVxEM9jXAWgSkSyAZgD2EAoRUQh8zrBJu4QDEXlKAaiUdwO6TlyeoHhbfvLhVHcD+BqAHSh03/Wp6s+KP/5y8dIH94iIvSokEVG5ys2w0+QXEG6G8UgWEVVyqP2gqi51/aB43ahbASwA0AvgByLyOwA+D2AvgHoAKwB8DsDfVLjJREQnlZlhp8wvINwM45EsIgq7u/DdALap6gFVzQL4EYBfU9UeLRgG8B0AyyLYEyJKonC7C0PLsNgeydqftd0l0+qOm9qNU14zNVeXw73Ddr6p5/sWmtrKlO1+6Z32kqnd0GS7TlpS9sii7eEJzjUf1r4B+4gzm8d351zRucv+3tFWUxvosY+1L2vn1nJ1l3ywc5Wp/c8d7za1D2/+oKl1NNium6vb3zC132nbbGqu59i1fb8atl08VzTY98Xfz7RdMi6vZ23n1QO99u9veccL4+67OoKaxZ5bsGvU1qbaqeQqFPpVkHcAuFJEmgEMArgBwCoR6VbVnuKlEW5Dwi99EGaGJSm/AH8yLOn5BTDDSj1QbAdZRBSiEDtzVHWliPwQwGoUJll+GYVD64+LyHQUOvbWAPhkeGslokTzNMM4yCJKugiulqyqXwDwhQnl60NdCRER4HWGcZBFRCc6boiI4snTDOMgi4gSM1UOEdUoTzOMgywi8jagiIgC8TTDYjvIuqn1VVMbUdtdMb/OdpysHbGdFI/1XGxqDelRU8s7DknuaZ1mapuz201tZtpun6ubJKj1A7NN7fiIbdc4nB4/91Vz3bBZZjRnt01abQdLd32vqbnmn3LND3Zu60FT++mmxfbxRm03zcACu1+XNNoOo2saTcm5fedn7Gu7P2eflxnpYHN1OecXc7xXPv76HePuz2uxc7q9Z9oGU7ukYbephebEhfxoUoWZYUnKL8CfDEt6fgHMsFJiO8giovCEOe8XEdFk8zXDzuhipCLynrA2hIiqKNyLkcYC84uohniaX2d6xfd7Q9kKIqqu8ucurAXML6Ja4Wl+lfy6UEQePdWPANhL9hJR7EiNHJ2aiPlFlAy+ZliQc7KuAfA7ACbObSDg3GNE8VdDXwE6ML+Iap3HGRZkkPUCgAFVfXbiD0RkU/ibdCqKnL7VZeGap8nF1ZkxpHa3bznLdvr0jdrOjNVH5praymN2PrDWlO1gPJqxXRj/3r/I1N7v6Dq6MGPbTrrr+0yto9nOm3V4YPx+1KfsvGQNjm6V7KB9nrYO2TnSMHWnKV3XOGRq26fY5S57m60912ufk/lNh0zN1Ynjsn7Edt1MTdn9Pebo7JqVtt+ou95Tro6droydc+2NLd3j7h+dbXfiAx1rTO2cOrttrnVWpqa+ApzIy/wCws2wJOUX4E+GJT2/AGZYKSUHWap602l+dm24m0NEVeHpp8AzxfwiSghPMyysYSRE5JdhPRYRTbIEdheOxfwiijlP8yvM62QFPPhJRN6pwYFTmZhfRHHmaYaFOcjydBeJ6LQ8vlryJGJ+EcWVxxkW2teFLiIyV0SeEZGNIrJeRD5drHeIyJMisrn4XzuvAxFNGtHybknA/CKKD1/zK8wjWa5h5CiAP1XV1SLSCuAlEXkSwMcBPKWqXxGRuwHcDeBzpR4+LacfE7o6FYbVzl3VnLK1z0zbYmq/tM0l6BmeamrrjtiOi53Hbe7u7bdzkA045up6x2V2Wy7K2HfGH03ban+3+Q1T+9mxS8bdz+Ztl8dQzr4V5rXbbqKc49PCwZztCOpyzJmVEdsR9KGWHaZ219Rgc1zZHhngqUH7fO4ftXOk7RjpstvSttqxDruWvrx9Y+yxu4bFjrnJuuceHnf/7dPt/t8y5ZippVD5HHGBJGTgdBpVzy+g8gxLUn4B/mRY0vMLYIaVUvaRLBFpK36S6xCRjjE/+tjEZVW1R1VXF/99DMBGALMB3Arg/uJi9wO4rewtJyIqE/OLiCZT4CNZIvIJAH8DYBBvjRkVwDkAoKrrSvz+fACXA1gJYKaq9hR/r0dEHBcuAURkOYDlADBvNueyJopKrX8FyPwiqm2+Zlg5f/l/BuAiVT1Y7kpEpAXAwwA+o6pHRYKdoKaqKwCsAIArLmv09CkkqgGenjQaIuYXUS3zNMPK+UhwbesAAB/2SURBVLpwCwD7xXUJIpJBIaD+VVV/VCzvE5Hu4s+7Aewv93GJKCTlXiMrnsMF5hdRrfI4v8o5kvV5AP8pIisBnLzOv6r+8al+QQof+e4FsFFV/2HMjx4FcCeArxT/+0g5G01EIYvnwKkczC+iWuZphpUzyPoXAE8DWAt3Y4TL1SicULpWRE5MZvQXKITTQyJyF4AdAD5cxnYAcHdITE3Z6wlmxHahXJgJNodSZ/q4qfUMtJnajt2208PJ8aXx1Gl2HXtH200tDzv3lcuyBrsfVzZsHHf/UN6uc3u77UwZytu3xyN9V5jaKyN2e29oGjG1323ba2o7Rm3n1HzH6+ji6pxxrfcB2+iCB7e+zdSG5mdMrbu+19QaxW7zs73nm9q2Yx2mtqRrfMfO3856ziyTQpOpufY1TL6ezxAir/ILCDfDkpRfgD8ZlvT8AphhpZQzyBpV1c+W8+Cq+gu4W6MB4IZyHouIIuRpQIWI+UVUyzzNsHLOyXpGRJaLSPcpWqCJKK48PqchJMwvolrmaX6VcyTro8X/fn5M7WQLNBHFU0Ku4s78IqpRPmdY4EGWqi6IckOIqIo8bX8OC/OLqMZ5mmElB1kicr2qPi0iH3T9fExbMxHFlaefAs8U84soITzNsCBHsq5FoSvnAyjshkz4b1VCytWFcyayOmpq59bZuZbqUrZDIr3fdnW4ZM7pN7W3z9ppalc22rmgto3aUfo/HninqV0yxXZ//H7bvnH3O1NTzDKddsosZzfIVY0vmZprjjTXfGtNYp/P2elgXSg/GWgxtaG8fd4/3NJnalc02udYHZ96frD5clPrarOvmcv+PjuvW9rxXrlwds+4+y2piOfzCkiibfypJi/zCwg3w5KUX4A/GZb0/AKYYaUEGWQdE5HPAliHt8IJ8HbcSERl8fh8hhAwv4hqnccZFqS7sAVAK4ArAPw/ALoBnAXgkwAWR7dpRDRpQu4uFJE/EZH1IrJORB4QkUYRWSAiK0Vks4g8KCKOYw+hY34RJUHI3YVhZVjJQZaqfklVvwSgC8DbVPXPVPVPUQitOcE2l4i8FuIgS0RmA/hjAEtV9WIAaQAfAfBVAPeo6kIARwDcFcGejMP8IkqIcD8khpZh5Vwnax6AsZejHQEwv4zfJyJPnWiBDnoLoA5Ak4jUAWgG0APgegA/LP78fgC3RbEvp8D8IqphIecXEFKGlXOdrO8BeFFEfozCWPD24kqIiE5S1d0i8jUUppwZBPAzAC8B6FU9eYb2LgCzJ3GzmF9EFEiYGVbOdbK+LCKPA7imWPo9VX25rC2P2EDezvvUnAp22kda7EE9V3fJwtYDprY247ieoWO4nM/bjpD6lO0K+tWQ/Rbjf+/6dVPb/OYsU2u5fNjUhlvHd+y8nrXrdMk42jUWZWw3zVXOJqmAHUti34Ku7pydI/bi3N/ZdpWpbZy7ztT+qmu9qd109gZT+8FaOx/Yrv5OU2tosc9xdk+zqY1k7HvgyELbFTWRa/9fzw6amuu1qFj5J412iciqMfdXqOoKABCRaQBuBbAAQC+AHwC4KZS1VigO+QVUnmFJyi/AnwxjfrnFIMNOmV9AuBlWzpEsqOpqAHYWTiKKr8o6cw6q6tJT/OzdALap6gEAEJEfAfg1AO0iUlf8JDgHwJ4Kt7gizC+iGlV+hp0uv4AQM6ycc7KIqFaF2124A8CVItIsIoLCZMobADwD4EPFZe4E8Eio+0BEyRVud2FoGcZBFhGFOshS1ZUonBy6GsBaFHJmBYDPAfisiLwBoBPAveHvCBElUoiDrDAzrKyvC4mo9gjCv5Cfqn4BwBcmlLcCWBbumogo6XzOMA6yiIjXPyeiePM0w2IzyBIAqTHfbrq6FxocnR45dczV5ejEcS339OBUU0s51rtwiZ2rq3/Ezuc0o9nOI7XhiO2weX7PAvt4m6aZGjpt99CRUdsl8qHNHxh3/zxHh9Hvd/7C1FydH7tGj5va9LTtfurL2207kEub2kX19nlKOb7Ffl/La6b2rdzVpvb8QUenlKM7Z9vxLlOTA3ZbGo46ZnZXu79tR+xigzMdm1J3bNz9h/vt6/qLo4tMbVnrFlO7IHPYrqASlZ34TmWYmF9AuBmWpPwC/MmwpOcXwAwrJTaDLCKKkKcBRUQUiKcZxkEWEXkbUEREgXiaYRxkEZG3h9qJiILwNcM4yCIibz8FEhEF4mmGRXqdLBH5tojsF5F1Y2pfFJHdIrKmeLs5ym0gohLKvUaWp2EWBWYYUQx4nF9RH8m6D8A/AvjuhPo9qvq1iNddliP5IVO7osHWrmw8ZGor6i43tQ393ab2m12rTC2r9iX40rrfMLX0sO0SyZkK8MohO19lOjW+o+htLW+aZfbmWk3thb7pprbLMQfXu1rtPFp7R23X0U8OXWpqrrnPPtr1gl1Hk+06+rPznzS1xw4uMbUPvG6fz9d2OuZN222f44xtpkHDUduhdWyu/bwy0mlfoed7zxt3v3fEdj81pu1z8lsdB+2GhPgZyddD7R64DzHIsCTlF+BPhiU9vwBmWCmRHslS1ecAhNSjSUSR8fiTYDUxw4hiwtP8qta0Op8SkVeLh+IdF1ApEJHlIrJKRFYdOOT63ENEYRAt70alM4z5RTR5fM2vagyyvgngXABLAPQA+PqpFlTVFaq6VFWXTu+0F7IkopDwSFY5AmUY84toEnmaX5M+yFLVfaqaU9U8gG+Bc5kRVRdPfC8LM4zIMx7n16QPskRk7BmVtwNYd6pliSh6UsEtyZhhRH7xOb8i7S4UkQcAXAegS0R2oTCj9XUisgSF8eR2AJ8I8lgK91xfpYy6elgcI9lpqcayH/uE5tSIqc1vsl08U1LDpnZDk23/OHDBc6a2aZ7t9tnQZyeXSjm+cB7Njx9Lf/nVm8wydWn7PM1otXOVzXTMX9abtZ0z63pt98uuQ/bUlaZG+9xd177J1PrzO+1yzbbDaEurfU6+u/YdpobDdv6unOMtMOw42+Zoxn420ZR93hv226+IVj954bj7Q7Pt/Gifu/pxU7uq8rdnMAk/OnUqYWVYpfkFBMuwJOUX4E+GJT2/AGZYKZEOslT1Dkf53ijXSUTl48nsbswwonjwNcN4xXci8vZTIBFRIJ5mGAdZRORtQBERBeJphnGQRZR0CkhlpwsREVWfxxnGQRYReXs+AxFREL5mWE0NstLiuCKF44l/cdgu5+p8yKkdGrvWMTPTF6i2tOGoXQnsij/ZvsPU+tpeN7VV09pMbcPQHFN7ofeccfe3bbFdM9khu1/bmm3HzYGuFlNra2o3tZ6tXaaGBvt83r7oFVNrTw+Y2r6c7RxaP2L3Y0bGPsfdXfa16HvZ/q6rO0cXHTe1c2fYObhe23KWqeWP2+6cTN/451ky9jl5eI+dS277kH0+396y1dQq5mlAJU2lGZak/AL8ybCk5xfADCulpgZZRFQZXz8FEhEF4WuGcZBFlHS8ijsRxZnHGcZBFhF5G1BERIF4mmEcZBElnMDfQ+1ERKX4nGEcZBGRt58CiYgC8TTDEjnIWubocnPNlZ0R+/S45h9rT9sOjs6UrU0NOL9YyrEtrt+9ocnOm3VZvZ2rtnnCnGOjF9nH7x22j//GNtvB0n9giqk1nDVqak6Ors7H31xsaq9MnW1qS6fZjqXtg52mtneg1dQO9NluogbHdHCOKdwwr+uwqX1q7tOm9n+n2G6aQyP2uTo81DTuft+E+wCwZdcMUxvI2rnKzm3cb2qVEvU0ocjJZlhy8gvwJ8OSnl8AM6yURA6yiGgMj08aJSIqyeMM4yCLiLw9n4GIKAhfM4yDLCLy9lMgEVEgnmYYB1lEFPqnQBE5H8CDY0rnAPhrAO0A/gDAgWL9L1T138JdOxElTZgZFmZ+cZBFRKF/ClTVTQCWAICIpAHsBvBjAL8H4B5V/Vq4aySiRAsxw8LMr9gMskaRx8HcW3NCdaXtnFRBuebvcnXEuDpxjuQHTe3XG2xbR3PKdlIE5VpvUK7n5ZLGnePuz5pl58I66pj46htDN5jaoa0dpta/xnbJNNqGIFcDFI4PTjW1N7obTG3bAbsOl5H9ttOlaa+dg6uh1/5FZqeIqR0dss/LBfV27q8rup80tfUjdt/SE6aKXzlwrllm24zppjbTMW/c/PoDplYRjfx8hhsAbFHVN0Xsc5wEE/MLCDfDkpRfgD8ZlvT8AhKRYWeUX47/7RFR4miZt/J8BMADY+5/SkReFZFvi8i0M9twIiJ4m18cZBEl3ImrJZdzA9AlIqvG3JY7H1ukHsAtAH5QLH0TwLkoHIrvAfD1yHeQiGpauRmGScyv2HxdSEQRKv9CfgdVdWmA5W4CsFpV9xVWU/gvAIjItwA8Vu6KiYiM8jJs0vKLR7KIqJIjWUHdgTGH2kWke8zPbgdgL/FNRFQmX/Mr0kFW8TvL/SKybkytQ0SeFJHNxf/ynAyiair3fKyAISUizQDeA+BHY8r/XUTWisirAN4F4E9C2YeIMMOIYsDj/Ir668L7APwjgO+Oqd0N4ClV/YqI3F28/7lSD1QHwbQx81/ltPIOFtfvphwdhy7NjvnA6mC7P3yyrGFiJ9KAWeaFYdv50XvUdvqkh213RSpr15l3NCdl+m1NO+07vskxp9nxo7ZLpnWq7ZQaabKvrabt6zM81e5Hzjb2YHDAFtcOzzS1W6YcM7WMo5umTcZ3Hl3d8JpZZkfbGlMbUvv+nJ6q/G9gIkfT0BlT1QEAnRNqHwt/TZG6DyFk2MT8AsLNsCTlF+BPhiU9v4DazbCw8ivSI1mq+hyAiTNU3grg/uK/7wdwW5TbQEQBRNtdGFvMMKKY8DS/qnHi+0xV7QEAVe0RETtdNxFNKl/n/fIUM4zIM75mmNfdhcW2yuUAMHe234e0iWJLUUl3IZXA/CKaJB5nWDW6C/edOEO/+N/9p1pQVVeo6lJVXTq9kyFFFJUIuwtrUaAMY34RTR5f86sag6xHAdxZ/PedAB6pwjYQ0Vg8J6sczDAi33iaX5F+XSgiDwC4DoWrq+4C8AUAXwHwkIjcBWAHgA8HfDTnnINjubp1XrKNHtgyYk+huKPVzufkXod9herOYDq2vvyQqU1N2U6UifOeAUDG8XxMS9lumolzibnmFtuZtXNrvW/hRlPb1m2Xa66zrTlbjtjlBtbYOcNcb/hzOg6Z2vy5E889Bs5r3mdqP9u/2NTWp2abmn3W4fyIM6vFPu8HRttM7WDObsuMdIupBZnXbUFd5XPaVeLE1ZLJCi/DSucXUHmGJSm/TlWrSoYxv5yYYW+JdJClqnec4kd21k4iqg5Vb89nqDZmGFEMeJxhXp/4TkSTI4rrZBERTRZfM4yDLCLy9lA7EVEQvmYYB1lESacA8p4mFBFRKR5nGAdZRMSOQSKKN08zrKYGWa7unbl1x01t80jlV65ocMz9tTtn56CanbZzRrm2z9WJ49Kayji2xdaCcHX6LGmwtd8860hFjw8AH97yPlN7c5ftLsxOsc/JzKvsPFr3nLXS1B7ut/Py7j421dRmzOwztQs7bDdNyvGl/sv755jaf/adZ2rXNL9h1xujSyP5eqg9aaLOsFrIL8CfDGN++cPXDKupQRYRVcjTzhwiokA8zTAOsojI20+BRERB+JphHGQRJR2v4k5EceZxhnGQRZRwhasle5pQREQl+JxhHGQREQLOlkFE5CdPM6zmB1nd6Smm9vbGNx1L2m4alx05O2tUR8q2Ybg6dubV2W1JOebods0PVYdgrR6u3x3W8fNyueYMOzdju4Rc2+ayO2e7aV4/ON3Uph527NegnTht94DtsHG5oH6vqX14/st2vWnbdbSwwf7uI0feZmrdrUdNbWaDrQ2pfX2O5O16+/K5cffn17WaZYLODxb09QnC10+BFDTDkpNfgD8ZlvT8AphhpdT8IIuISvD4fAYiopI8zjAOsogSz9/JVYmISvM3wzjIIiJv25+JiILwNcM4yCIibz8FEhEF4mmGcZBFlHQKOGbkICKKB48zLJGDrEWZYJ04ObWv2nRHJ06TYw6uFkczjevxYJtT3B0XjuVcj+eaX2yioPONBe0Q2Zq13SXHd7eYWn2n3ba6QfvpY/Ne29XTf96wqV1U32Brna+dcjtLefKo7WI6nq03tYzYDpvW1Kipbcra57k33zyhcsgscyBnf68jbfe/Mczj455+CiS3IBlWq/kF+JNhzC+AGXZ6iRxkEdEEfuYTEVEwnmYYB1lE5O01ZoiIgvA1wzjIIiJvD7UTEQXiaYZxkEWUdApvp6QgIirJ4wyr2iBLRLYDOAYgB2BUVZdWa1uIkkyg3h5q9xkzjMgPPmdYtY9kvUtVD4b1YFuydq6l5pR94menK59rqTdvuzBa6myXiIurmyZox45L0E6cBkf3UKV+aac+w+NHLzW12QsPmFr/PNvpknO8Pu+dtd3UVvRebGoLGvab2u1T+uwGOuwaPW5qFzTtMbW5s23nzEDe7scXd/+Gqb3cM9vUcvnxr9nIgH1tGpptl9C0FvvePm/qqf507j1F/TRCDigROR/Ag2NK5wD4awDfLdbnA9gO4LdU9UioK59cXmQY8yu4MDMs6fkFRJFhFeQXEGqGhZlf4c3OSETxpVrereTD6SZVXaKqSwBcAWAAwI8B3A3gKVVdCOCp4n0iojPjaX5Vc5ClAH4mIi+JyPIqbgdRsp04n6GcW3luALBFVd8EcCuA+4v1+wHcdmYbX1XMMCIflJth5Tmj/Krm14VXq+oeEZkB4EkReU1Vnxu7QDG4lgPAvNnV/maTqHZFfD7DRwA8UPz3TFXtAQBV7Sn+/cfVaTOM+UU0eSLMsDPKr6odyVLVPcX/7kfhMNwyxzIrVHWpqi7t6nRcgpiIwlH+14VdIrJqzM15JEdE6gHcAuAHk7k7k6FUhjG/iCaRp/lVlY9XIjIFQEpVjxX//V4Af1ONbSGiYOcpTHAwYDfdTQBWq+q+4v19ItJd/BTYDcCe/RsDzDAin5SdYZOWX9U6hj0TwI9F5MQ2/H+q+tMzfdDjanfnuJ2mCd3pYF/KurpfNmQ7TW1OnaNdJaCe3KCp7XPM+3Rhxu5Ic8p2ibhMnEvM1YkUdB6xZQ12uefSdv+Xz/8PU5ufsZ0kHWnbcTI9ZdexJ2df2zVD80xtf67H1Gal7RxkU1J2365vfsPUpjrment20M5N9i8vXWtqHc/brq3M8fFBoI4DHP1n2dd/b/cUU/ujm56xvwzg/zirp6EAcpEdar8Dbx1qB4BHAdwJ4CvF/z4S1Yoj5n2G1Wp+Af5kWNLzCwg/w8rOLyDKDDvj/KrKIEtVtwK4rBrrJiIrivMZRKQZwHsAfGJM+SsAHhKRuwDsAPDh0Fc8CZhhRH4JO8PCyi+ejUlEkUxJoaoDADon1A6h0K1DRBSekDMsrPziIIso6RRA3s+rJRMRleRxhnGQRZR4FZ34TkTkCX8zjIMsIvI2oIiIAvE0w2I7yOrL226QRrFPcodjbqmco4Mn6Dxa722y63XO3xVQe8puy2P9C0xtenqjqc1zdOdsG7XdLhOfl+607fJw7X9/ftjU/m3gLFPrzvSa2m+3BuvMP5izz92A44/l0nq7r7PSW00tI3bitKza+dqmpZocNbt9Lw7b7Tucc3T7tNn3Rf0xO6fXtOd3jbs/ssB1LTvb1TMw327Hzc12rrKKeRpQtSzMDEtSfgH+ZFjS8wtghpUS20EWEYXE4/MZiIhK8jjDOMgiSjwFzuBoBhFRdfmbYRxkEZG3h9qJiALxNMM4yCJKOo8PtRMRleRxhnGQRUTefgokIgrE0wyL0SBLx3XBpGG7MBrFzo/V5ehCcdmYDTZ/10WZZlu0m4L9Odsl05my8zm1pGwXxu+22TmoMmJ/12WKY1s6HZ0oE7k6jBrEvj1+c8oBU3th2C6XVftaNIjtVulK2+dz/YjtCHJxTKOFYzk7kdaS+mBvc9d8aP/n0K+Z2hOvLza10X67b62OTqHRHeO7czJt9v3Z1Gz3ITVg92HA8RxXzNOAqh1q/sbCzLCk5xdQnQxLen4BzLBSYjTIIqJo+HshPyKi0vzNMA6yiJJOAeT97MwhIirJ4wzjIIuIvP0USEQUiKcZxkEWEXkbUEREgXiaYRxkESWeetv+TERUmr8ZFqNBloybm6pFbFdLk6N7LahFdXZuKddcWK55pFzLzXB0nQSdI6zZMadXUEFmMAu6HRmxHSIuZ9cdM7Vhx/yQQR9vejpYx8mCOvscv54dDPS7rtdxFHa9adg/3PxB+95r226f+Xyd/d3hDywbd19TtoOnd4F9nlrOPWJqW7KOTrFKKKCeXi25dojJiTAzLEn5Vc62VCPDkpRfADOslBgNsogoMp5+CiQiCsTTDOMgi4i8PZ+BiCgQTzOMgyyipFP1tv2ZiKgkjzOMgywi8vZTIBFRIJ5mWNDzDEMnIjeKyCYReUNE7q7WdhARoPl8WTdihhH5xNf8qsqRLBFJA/h/AbwHwC4AvxKRR1V1w5k8rqtLxtWF4lrOVQu6DpeUa/zqmJcraJeMa72udbjm0jJzpgXch6Dm1Nm5qw465j5rCvhauDqbXFz7f0HGbotrTq9/6j3P1FrTtrNnMG+7vbTFdvY0HbTLjTbZF3z3jdnxj7XPMafbdLsdf3fhE6Z2SX2wOdJK83dKCl/5lmFJyq9TPd6ZqDTDkp5fADOslGodyVoG4A1V3aqqIwC+D+DWKm0LUbIpCp055dyIGUbki3IzbBJV65ys2QB2jrm/C8A7qrQtROTpNWY8xgwj8omnGVatQZbjwLO9apqILAewHADmzeY5+kRRUFVoLtgFFOmkkhnG/CKaHD5nWLW+LtwFYO6Y+3MA7Jm4kKquUNWlqrq0qzPYlXuJqHya17JuVDrDmF9Ek8fX/KrWx6tfAVgoIgsA7AbwEQAfrdK2EJGnh9o9xgwj8omnGSZapTPyReRmAN8AkAbwbVX9conlDwB4E0AXgIPRb2HkamE/amEfgNrYj7H7cLaqTg/6iyLy0+Lvl+Ogqt5Y5u/UlHIybEx+AbX3fouzWtiPWtgH4K39KCu/gIoybNLyq2qDrEqJyCpVXVrt7ThTtbAftbAPQG3sRy3sQ1LUwmtVC/sA1MZ+1MI+ALWzHxNV7WKkRERERLWMgywiIiKiCMRxkLWi2hsQklrYj1rYB6A29qMW9iEpauG1qoV9AGpjP2phH4Da2Y9xYndOFhEREVEcxPFIFhEREZH3YjXIiuus9yLybRHZLyLrxtQ6RORJEdlc/O+0am5jKSIyV0SeEZGNIrJeRD5drMdmP0SkUUReFJFXivvwpWJ9gYisLO7DgyJSX+1tLUVE0iLysog8Vrwfu31IGuZX9dRCfgHMsDiKzSBrzKz3NwFYDOAOEVlc3a0K7D4AE6/JcTeAp1R1IYCnivd9NgrgT1X1QgBXAvjD4vMfp/0YBnC9ql4GYAmAG0XkSgBfBXBPcR+OALiritsY1KcBbBxzP477kBjMr6qrhfwCmGGxE5tBFmI8672qPgfg8ITyrQDuL/77fgC3TepGlUlVe1R1dfHfx1D445iNGO2HFvQX72aKNwVwPYAfFute7wMAiMgcAL8B4H8X7wtitg8JxPyqolrIL4AZFkdxGmS5Zr2fXaVtCcNMVe0BCgEAYEaVtycwEZkP4HIAKxGz/Sgeol4DYD+AJwFsAdCrqqPFReLwvvoGgD8HcGIeiU7Ebx+ShvnliTjnF8AMi5s4DbJKznpP0RORFgAPA/iMqh6t9vaUS1VzqroEhQl9lwG40LXY5G5VcCLyfgD7VfWlsWXHot7uQ0LxNfJA3PMLYIbFTbUmiK5EyVnvY2afiHSrao+IdKPwqcRrIpJBIaD+VVV/VCzHbj8AQFV7ReTnKJyf0S4idcVPUb6/r64GcEtx3rxGAG0ofCqM0z4kEfOrymopvwBmWFzE6UjWyVnvi10HHwHwaJW36Uw8CuDO4r/vBPBIFbelpOJ35vcC2Kiq/zDmR7HZDxGZLiLtxX83AXg3CudmPAPgQ8XFvN4HVf28qs5R1fko/A08raq/jRjtQ0Ixv6qoFvILYIbFkqrG5gbgZgCvo/Ad9F9We3vK2O4HAPQAyKLwifYuFL6DfgrA5uJ/O6q9nSX24ddROHz7KoA1xdvNcdoPAJcCeLm4D+sA/HWxfg6AFwG8AeAHABqqva0B9+c6AI/FeR+SdGN+VXUfYp9fxf1ghsXsxiu+ExEREUUgTl8XEhEREcUGB1lEREREEeAgi4iIiCgCHGQRERERRYCDLCIiIqIIcJBFREREFIE4XfGdJomIfBFAPwpX4n1OVf/9DB5rCYBvFh8rB+DLqvpgGNtJRDQR84t8wkEWnZKq/nUIDzMA4HdVdbOInAXgJRF5QlV7Q3hsIiIn5hf5gF8XEgBARP5SRDaJyL8DOL9Yu09EPlT893YR+XsR+aWIrBKRt4nIEyKyRUQ+earHVdXXVXVz8d97UJgbbPok7BIRJQTzi3zFI1kEEbkChTmkLkfhPbEawEuORXeq6lUicg+A+1CY6LMRwHoA/xxgPcsA1KMwrQgR0RljfpHPOMgiALgGwI9VdQAARORUE9eeqK8F0KKqxwAcE5EhEWk/3SH04gz33wNwp6rmQ9x2Iko25hd5i18X0glBJrEcLv43P+bfJ+6fcsAuIm0AfgLgr1T1hYq3kIjIjflFXuIgiwDgOQC3i0iTiLQC+EBYDywi9QB+DOC7qvqDsB6XiKiI+UXe4teFBFVdLSIPAlgD4E0A/xHiw/8WgGsBdIrIx4u1j6vqmhDXQUQJxfwin4lqkKOsRERERFQOfl1IREREFAF+XUihEJFLUOi+GWtYVd9Rje0hIgqK+UVR4deFRERERBHg14VEREREEeAgi4iIiCgCHGQRERERRYCDLCIiIqIIcJBFREREFIH/H6PNfUOHNnjcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=[10, 10])\n",
    "ax1 = plt.subplot(221)\n",
    "ax2 = plt.subplot(222)\n",
    "xds_c_read[0,:,:].plot(ax=ax1)\n",
    "xds_f_read[0,:,:].plot(ax=ax2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f2eaddc8510>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAEmCAYAAABVkVXzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5gcVZk/8O87k8l9cp0EhsxAQggXUbk4xADKNSKgCLqyi64uKmt0V1bw8lNYXMVd3dVdFS+gbhSWuA8risIDcicSZFm5mISAiQGSACFDQkJIJvfMZGbe3x/dwZl5T89UdVd1ndP1/TxPP6TfdHed6pl8OV1dbx1RVRARERFROuqyHgARERFRLeNki4iIiChFnGwRERERpYiTLSIiIqIUcbJFRERElCJOtoiIiIhSVJXJlojUi8iTInJn8f4MEXlcRFaJyC9EZHg1xkFEFBfzi4gqJdW4zpaIfBZAG4BxqvpuEfklgFtV9WYR+TGAp1T1R4O9Rv2ocdrQOPX1+0e3jDOPWdG+3dT4OPs4n8eWt8clvU0A6Hx1zWZVneL8S4dWGaV70Rv14QCAzei6T1XPjvWkQCWRXyMaJ+jYpuZ+tekTRyU2xo7OHlMbViemNrYh2c/Xu/bZ35vVL71qH6j2cccc1mwfF9HeHvv/rb3ddhv2HXDr7rWvN3nUsCGft25bp6lt277X1Orq7ft+1EGNpuYar+tfputxUf9P7voNeOr5TaZWP8x+hnjjwRMibcM15q17uk1tU8ceU9vzyupY+QXEz7As8iv1yZaItABYAODrAD4L4DwArwI4UFW7ReREAFer6jsHe52RUw/TQ/7qW6/ff+rf7cOP+cJ9psbH2cf5PLa8PS7pbQLAc9e+d4mqtjn/0mGqjNALJd7/+H6oa2NtI1RJ5dfkGUfpO7/6s361Gy58Y2LjvGv1VlObOKrB1E6aNjaxbQLAHzbsNrULLv2xqfV02f+pvvKbq8re7pqOfaa24tWdptbgmHC6vLa7y9Q+9Kapjkf295m7njO1e+59xtQaJ9n3/fdXn2FqDY7hdjrmEK7dcswXnUY4Zlstf/k9Uxsz5WBTe/a690bahmsy/Ks/2Un4D25dYWpPf/NdsbMlboZlkV9DT90r910AXwCwfxo/GUCHqu6f5rYDmFaFcRCRgwCol6jHAIrys/AE84vIc7EzLIP8SvWcLRF5N4BNqrqkb9nxUOeui8g8EVksIot79ri/MiGiytVLvFseJJlfe3d0pDJGIirwPb/SPrJ1MoD3iMi5AEYCGIfCJ8UJIjKs+OmwBcB615NVdT6A+UDha8SUx0qUS2Ud2cqHxPJr8oyjmF9EKQkhw1I9sqWqV6pqi6pOB3ARgAdV9a8BLALw/uLDLgZwe5rjIKJBxDyqlZcjW8wvokAEkF/VOGfL5YsAbhaRrwF4EsD1Qz3h6JZx+L8+Jwn7coJziI/zeWx5e1zS2wSAUdc6yyWF8KnQM7HzSwH0DDiD+WO3LI+0sQuPbxnyMXWOn1/Uk+HvWZPsV5w//fePmVpvwo1Ya7baE/Ojngzv8uAzthvviRdt08HAn6HLmXMPN7V3v/FAU3OdDO/yoZuWmdo6x9ge+6fTo72gw471a0xtz9aNjkdGO0F+pGNG42o4cJxHj09+M9Im+gkhw6o22VLVhwA8VPzz8wBmV2vbRFRaIaiyHoXfmF9E/gohw7I6skVE3hDvPxUSEZXmf4ZxskWUcwKu20VE4QohwzjZIiLvPxUSEQ3G9wzjZIso5yRHHYZEVHtCyLCqrI2YhIHLXUTpCgGAi9paIz2u3tHJctaM8aaWdOdOVI0j7Lz4bS2228i1tMO2Aeumrd5il87YtteuW1UJ1/t5y9J2Uxs+rD7Sc79/3hGRtutz56HrcT4s19NSP1IvG2WX5hjMF3atysVyPUlxLdeTZIblKb8AfzKsGvk1coxdo/BhR+ehaxmeHzzxsqldu+AJU9vw5EJTqx9u1+7c+tC3TK0S3Y5/Ao1jRsfOlrgZlkV+8cgWUc6F0MlDRFRKCBnGyRZRzoVwjRoiolJCyDBOtojI+0+FRESD8T3DONkiyrnCyaWeJxURUQkhZBgnW0Tk/adCIqLB+J5hnGzF5Oo0idpVVImxw23Xi4urI2XJ+h2JjSNq15NL1MdF9fFfrzC1vbvsmmnv/OHjplZqXcGBRjXaTqDD/+5XpvbKU4tMbcIp95ra6vuG7uaJOrZk10aM9xwKU97zC/Anw6Lm14HTJ5ia631yuem+50xNe213p4vrcZv32FrTqGg/W5dhCeVOCBnGyRZRzkkKS12IyEgADwMYgULO/EpVv5LoRoiIkHyGpZFfnGwRURqfCjsBnKGqO0WkAcAjInKPqj6W+JaIKPcSzrDE84uTLaKcS+PkUi1cLXln8W5D8RbGFZSJKChJZ1ga+eX72o1ElLL95zvEuQFoEpHFfW7zzOuK1IvIMgCbADygqvbkOSKiCsXNMGSQXzyyRZRzAmC446ThIWwearkLVe0BcKyITABwm4i8UVWXlzlMIiKnMjKs6vkVzGRLkW7XTCXdKNVYb+zYqXadqiy8pbmx7Oeq4zDvUxtt901nt10g7fblr5jao79/yW7D0UFz5pubow7Rbvdzp5jaqV+8y9QOPMauVSZ1tktn24DF3yrp5ElSmteoUdUOEXkIwNkAcjnZYn75kV9AshmW9/wCaj/Dksovfo1IlHMiQJ1IrNvQrylTip8IISKjAMwF8EzKu0JEORQ3w4Z+veTzK5gjW0SUFoEk347YDGCBiNSj8KHul6p6Z9IbISJKIcMSzy9OtojyToC6hCdbqvo0gOMSfVEiIpeEMyyN/OJkiyjnBIDU84wCIgpTCBmW6mSr1FVYReRGAKcC2FZ86EdUdVmaYyGiEgRpfI0YPOYXUSACyLC0j2w5r8Ja/Lv/p6p2kbkSJo0ejovaWhMZ1Mhh0WbASx2dJlv2dCcyhsEcMNaux+eybNMeU3N1/TQMmPHv67EdJVEt37TL1PY5uqxc2xg4jlKPu9PRueMye479ffi/h1aZ2mdOPjjS6331wRdM7eK2FlMbO2GMqTWMGGdq379ktqkdPL6h3/1uR4Naj6O28PmtjtdKqMNLJPGvEWuEl/kFRMuwWs0vwJ8My3t+AcywqFKdbPEq0kRhkDq/D8FngflFFA7fMyz10Q1yFdavi8jTInKNiIxIexxE5CbFk0vj3PKC+UXkv7gZloXUJ1uq2qOqxwJoATBbRN4I4EoARwI4AcAkAF90PVdE5u2/nP72ra+lPVSi3JJ6iXXLC+YXURh8z6+qHXdT1Q4ADwE4W1U3aEEngP8CYL8cLjxnvqq2qWrbuImTqzVUonwRgdTXxbrlDfOLyGMxMywLqW611FVYRaS5WBMAFyCnS3gQ+UDArxFdmF9EYYibYVlIuxvReRVWEXlQRKag8B4tA/DJlMfRz2GTRpvaorXbTW2vY42rpLk6bY4/wI7PpaUxWtfP3On9O03uf2GbeYxr3bbDJttxrH7NdjhFFbWD6N1vPDDS425evM7U5p51lKl97s5nTW24433vcozvO7973tTEseDpo186rdQw+1m7Y1+/+x/56R/MY1550Xbt9HZ3mZprHbWyiHufyM/8AmyG5Sm/AL8zLE/5BTDDokq7G9F5FVZVPSPN7RJRHIK6HH41OBTmF1Eo/M8wXkGeKO8CuCAgEVFJAWQYJ1tEOScBBBURUSkhZBgnW0Tk/SF4IqLB+J5hnGwR5Z3k69pZRFRjAsiwXE62etR2rVSjc8cl6jqNLht22s6QZRvtul/jR/T/Mbu6dkI0dqRdp6ve0ZHi2t89ji6Yay+wnUBX3bfa1D563pGm9syWTlM7cpK9sPghjf3HfOzhTeYxizbuNLV9nXa/ho+y+18OAVDneScP9Tcww/KUX0BtZFgt5BfADIsql5MtIupDkMsLlRJRjQggwzjZIqLcXKiUiGqT7xnGyRZR3gVwvgMRUUkBZBgnW0Q5JwEcgiciKiWEDONki4i8PwRPRDQY3zOMky2ivAtgXTEiopICyLBcTrZWvmrbi7MyooJDn1v22NZp10Kpu/cNvdina8HWSpwzc0Kir7fwRbvQ7p6ublM7ZPIYU/vymYea2mV3PGNqt63cbGr/etZMU1u2aY+pXfLjx0ztt188xdRGDvj09Y1zDjePGfEuW9vbY9u/RdzhMuG7znJJEsC6YtSfLxnmS34BfmdYreYXkHyGxc0vIIwMy+Vki4j6qAPqhtdnPQoiovIEkGGcbBHlnkDq/P5USERUmv8ZxskWUc6JJL+umIi0AvgZgAMB9AKYr6rfS3QjRERIPsPSyC9OtohyT9Jom+4G8DlVXSoijQCWiMgDqvqnpDdERHmXeIYlnl+cbBHlXQrXqFHVDQA2FP+8Q0RWApgGgJMtIkpWwhmWRn4FM9mqrxOMH/nn4brWIW1wtH66Ol4qMboh2kl4UTtoDp04MtLjNu62r3fQOLtQqC/uf2GbqbkWVHUtvHrWjPGmNnf6OMdWDjaVO/64wdTEsfD49887wtS+eK9dtPV9R042tTdPtV1Pt1x2sqkNd/7j7z+WERHzwdUFlJx0z3cQkekAjgPweGob8dzA/AKyyTDmV3RRMizv+QXUfoYllV/BTLaIKCUCSH3sTp4mEVnc5/58VZ1vXlpkLIBfA7hcVW3/OxFRpeJnWNXzi5MtopyT8s532KyqbYO+rkgDCkF1k6reWu74iIgGU0aGVT2/ONkiyjsB6hI+BC+FqxVeD2Clqn4n0RcnIuor4QxLI79SvTCFiIwUkSdE5CkRWSEiXy3WZ4jI4yKySkR+ISLD0xwHEQ1O6uti3SI4GcCHAZwhIsuKt3PT3YtkMb+IwuF7fqV9ZKsTwBmqurN4SO4REbkHwGcBXKOqN4vIjwFcAuBHKY+FiFwk+Us/qOojAPxerGxozC+iECScYWnkV6qTLVVVADuLdxuKNwVwBoAPFusLAFyNIcJqbEMdTpo2dtDt/fHVvaYWtZMnapfOp6/9faTHff/Sk0zN1eGzfJNd42y6o8Ony7GulM9cnYcNjn8MrrXQXF1AUTt8XLXP3PWcqV3jWLvrm2cfZmoudY6OnIPGuP4pDf0z643477nX0ZFUV2JtxLgE8P7qy1modn4B5WcY8yt5UTIs7/kFMMOiSv2cLRGpB7AEwGEArgOwBkCHqu5fhbMdhetXEFEWUjiyVSuYX0QBCCDDUp9sqWoPgGNFZAKA2wAc5XqY67kiMg/APABobW1NbYxEuZbCRU1rBfOLKAABZFjVRqeqHQAeAjAHwAQR2T/RawGwvsRz5qtqm6q2TWlqqs5AiXKorr4u1i1vmF9EfvM9v9LuRpxS/EQIERkFYC6AlQAWAXh/8WEXA7g9zXEQUWkihasvx7nlAfOLKAxxMywLaX+N2AxgQfG8hzoAv1TVO0XkTwBuFpGvAXgShetZEFFGfD8EnxHmF1EgfM+wtLsRn0ZhTaGB9ecBzK7otR1dDBt2dpb9eq4undEVrN0VtevH5Y5/OqPs56btsMl2Xa3Vr+2O9NyRw+w/hgPG2EsUtW+3HVn3rOkwtXNmToi03W27u0ztu4+1m9rlc1oivV5Urt/RgeucvdBhx7Zg8TpTe22XfZxrXbayBHByaRbSzC8g2QxjfkWXZIblPb8AZlhUvII8EeXmq0Eiqk2+ZxgnW0Q5JyKoi78QNRGRF0LIME62iMj7Q/BERIPxPcM42SLKuwCuUUNEVFIAGcbJFlHuiffnOxARleZ/hgU72XJ1RbjWsnK55PM/MbWpRxxjalf+la196p9uMrVJB88ytVdXPR1pLC7v+Rdbc3X4ZLHeWNSunfEj7a9WlLXhAOBNU+zaaj9cbK8b2a22m2eYo7ll1HA7lkvfmn7njuNX1KwiNmOC7WaaNmmUqe3o7Da1qL/vQ5EAPhXWoiQzjPkVXdoZlqf8AphhUQU72SKihATQNk1EVFIAGcbJFhF5fwieiGgwvmcYJ1tEeScCqfO7bZqIqKQAMoyTLSICPA8qIqJBeZ5hnGwR5Z0IpKEh61EQEZUngAwLZrK1vasHi9Zuf/3+voS6GPbbs2OXqbU41hZ7wyknmNrqxc+ZWv0I243Ru8+uDaW9Pabm6gQ68W+idQf97r8ui/S4KKJ27bjs3teb2DgA4O/bDir7udedf2SCIwE27LJdNV9buNrUhjvWgrzmXYf3u18H+3t8yXHNpvblzfZnMXp4Up/kxPtPhaEbmF9AshnG/HLzJcNqNb8AZlhUwUy2iCglAu+DioiopAAyjJMtopyTAC4ISERUSggZxskWUd4F8KmQiKikADKMky2i3PP/fAciotL8zzBOtojI+0PwRESD8T3DgplsKfp377jWVOp1LebkMHnmm0xt1+ZXTO0NTbYj57Z5s03t0iljTG3Rbx43tRFTDja1HettF0hUrk6gUz/6vbJe6/pvfbzscbjs67GdPAtf3G5qbQfZtcbWbbNdT3eu3Ghqcw6ZaGq/WmbXIDvziCmm9r4jJ5tad8TmsAPGRGsxbhwx9D8v1zb/43/XmtrksXb9sc+cZH+fAODrQw+tP/H/U2HoBuYXkGyG5T2/gGwyLO/5VWq7lWRY7PwCgsgwv6eCRFQFxaCKcxvqFUVuEJFNIrK8CjtARLkWM8OGerUU8ouTLaK8E0Dq62PdIrgRwNnpDpyICLEzLIIbkXB+BfM1IhGlRYCEz3dQ1YdFZHqiL0pE5JRshqWRX5xsEeVdAG3TREQlBZBhqX6NKCKtIrJIRFaKyAoRuaxYv1pEXhaRZcXbuWmOg4gGI5C6+lg3AE0isrjPbV7We5E05hdRKOJlGDLIr7SPbHUD+JyqLhWRRgBLROSB4t9do6rfivpC44YPw9wZE16/f9+arWUP6rU1fzQ1ccyK73Vs44LDJ5naTy+03Y1fGm87gR5fucmOZcxIUxvjWNNs04u222jv9ldNradzj6mlbXSDfe9277OdRq7unu88/KKpbdrRaWr1dWJqv9pq93Xn3n2m9lJHtPfkwRe3mdotS9sjPXfsSNvh89V3zLQPHNBtNszuFtod+9Uy0f4+udYkK1v8Q/CbVbUtuQF4KbX8ApLNMOZXZcrNsLznFxBshlU9v1KdbKnqBgAbin/eISIrAUxLc5tEFJOI83/Wecf8IgpEABlWtW7E4slmxwHYfwGXS0Xk6WKLpb3gSOE58/Yf5tu82X4KIqIkpHLph58DeBTAESLSLiKXpL4bKWJ+Efks8Us/JJ5fVZlsichYAL8GcLmqbgfwIwAzARyLwifHb7uep6rzVbVNVduamuyF3YgoAYLCIfg4tyGo6gdUtVlVG1S1RVWvT39H0sH8IvJc3AwbQhr5lXo3oog0oBBUN6nqrQCgqhv7/P1PANyZ9jiIyE0gUa89kzvMLyL/hZBhqU62REQAXA9gpap+p0+9uXg+BAC8FwCvMk2UlQDaprPA/CIKRAAZlvaRrZMBfBjAH0VkWbH2jwA+ICLHorBk2IsAPpHExurE0RbhMHqyPcfV1Rnzma/dZmrXth1tags/PcfU/vUs28nx4hy7DtTqLbZrY9Zk27VxwGj7o5pz1QOm9uqqoXPftSZZJU49uNHU7lnTEem5nzzRvicHjYn2a/nt368ztbWv7Ta1y+e0RHq9s2aMN7WbF9ttuNa027bbrod27eO2E+gfZg99frVrLbQ3TLVrsPUi2u/70PxfVywjVc0voPwMy1N+Af5kGPPLjRlmpd2N+AjgfDfvTnO7RBSPJHwF+VrA/CIKh+8ZxivIE+Wd+P+pkIiopAAyjJMtIgLE70+FRESD8jzDONkiyj3xPqiIiErzP8M42SIiqOdBRUQ0GN8zLJzJlhS+lo2rafRwU1v4nQtN7W1/+0NT69yxxdRW/u73pnZMu12T6ql/f6epTR9n159y1body0W51p/att52mkTp0nEta3DJ539iatd/6+NDvlYprrXAXF0wT72y09S+5FjPa9Rw+6v62VNtx9Qh4+376fLlhc+b2j/PPdTUbrjwjZFe7wdPvGxqT661a9PdNWl0v/vnzrJr1R3YaNeWu+sZuy7dESfbTqiyCLz/VBi8MvMLiJZhecovwJ8My3t+AcywqMKZbBFRSiqYCRARZc7/DONki4giLWFBROQtzzOMky2inFMRaB2jgIjCFEKG+T06IqoOz893ICIalOcZxskWUe753zZNRFSa/xkWzmRLtV8nSOMIO/R9Pb2mdkKz7Z54pN12kPzk3y42NVd3i6tbZvsrtvvkH+9fY2p/+1bbeXGoo/vE1bnjsm+v3Y8okl5XzMW1xlsPHG1KDhceb9cCc63xNf+xtab2zbMPi7SNj5zQamqL1m43tdMPGRfp9T412475bx3dPL9c0v935e4VG81jXF1Qsw6w64qJRns/o/C9bTp4A/ILSDbD8pRfgN8Zlqf8AphhUYUz2SKi9HgeVEREg/I8wyoanYi8I6mBEFFGROLfagDzi6hGBJBflR7Zuh5AQlclI6LMeP6pMCXML6Ja4XmGDTnZEpE7Sv0VgMnJDoeIsuD7+Q7lYn4R5YPvGRblyNbbAXwIwMCzGQXA7MRHRERVJt5fELACzC+imud/hkWZbD0GYLeq/m7gX4jIs8kPKZqTWhpNrS5it5ur48fFtf6Wqwuma5ddW+yWXy02tcdX2rWhnnt0iam94y9OMzXXGld1w+yaab37ukwtC3On2y6Ye9Z0lP16F7XZ7ptzZk6I9Nxv/952As2aYjtj1m7dbWpRu3mi/u4NtHPvPlMb71gL70PHNJf1+pEEsK5YBbzMLyD9DGN+VSbJDKvV/AKYYVENOdlS1XMG+btTkh0OEVWf/9eoKRfziygP/M+wxEYnIo8m9VpEVGVSF+9WY5hfRIHzPL+SvM7WyARfi4iqyPeTS6uA+UUUMN8zLMnJVnKXgiWi6hH/D8FXAfOLKFQBZFiqoxORVhFZJCIrRWSFiFxWrE8SkQdEZFXxvxPTHAcRDSHhiwKKyNki8qyIrBaRK6qwB4ljfhEFJOGLmiadYUke2XLtQTeAz6nqUhFpBLBERB4A8BEAv1XVbxR34goAX0xwLK/r6LJdO09tsGtIfXqO7Rb5wCXnmdrPr/+Nqbk6fHa9+pKprXzI1lyWLllvi45unpf+5+9M7cSvPWRq655+ashtutZWi+q5rbaD6PCJthvF5dgDbVdN85jyfy0/f/cqU+tydG6teNl2YH3i5BmRtrGmw3bfbNzVaWquDqSB6+PdstSuNfb9846INI7kJPupUETqAVwH4B0A2gH8QUTuUNU/JbaR5HmZX0C0DMt7fgHZZFje8wtghkUVe3QiMq74yW6SiEzq81cfHvhYVd2gqkuLf94BYCWAaQDOB7Cg+LAFAC6IPXIiSoxKXazbEGYDWK2qz6tqF4CbUfg3nznmF1FtSjC/gBQyLPIUXEQ+AeCfAezBn89vUACHAoCqLh/i+dMBHAfgcQAHqOqG4vM2iMjUEs+ZB2AeALS22lk2ESUk/qfCJhHpezGm+ao6v/jnaQD6XhyoHcBbKxhdxZhfRDUuXoYNll9AChkW53jn5wEcraqb425ERMYC+DWAy1V1u0T8zrS48/MB4Pjjj+cJrEQpUBFo/MVZN6tqW4m/c71Y1v9+mV9ENaqMDBssv4AUMizOVHANAHuJ2iGISAMKQXWTqt5aLG8Ukebi3zcDsJcmJqLqUEBj3obQDqDvoZwWAI4TeaqK+UVUq5LNLyCFDItzZOtKAL8XkccBvH42nap+utQTpPAR8HoAK1X1O33+6g4AFwP4RvG/t8cZNBElSdEbMYEi+gOAWSIyA8DLAC4C8MEkN1AG5hdRzfI/w+JMtv4TwIMA/ggg2uKCwMkonHj6RxFZVqz9Iwoh9UsRuQTASwAujDEOAMDKzXtN7eimEaY2tsEevLt8TovjFe0P6pCm0aY29y/fEW2AFRg+rPyuike/dJqj2r+2qsN23qzfbrtRXL+8e7vtj/6lbfZn4erkca0FtnSjPdhQSTfPt86dZWof//UKU7vwePs7sGWP7dKpZD3HKL569pFlP7fXeaS7PEnGlKp2i8ilAO4DUA/gBlW1P4Tq8iq/gKQzLD/5BfiTYXnPL4AZFlWc34puVf1snBdX1Ufg/u4TAM6M81pElA4F4Ojoruw1Ve8GcHeyr1oR5hdRjQohw+JMthYVu2t+g/6H4bckNRgiyoYmewjeR8wvohrme4bFmWzt/77yyj6111uniShMaXwq9BDzi6hGhZBhkSdbqhrt8rREFBzPc6pizC+i2uZ7hg052RKRM1T1QRF5n+vv+7RDE1GI1P9PheVifhHlQAAZFuXI1ikodPGch8LkUQb8N5OwcnXtuAyL2OzQ7fhBzWubZmqLX0j/FI/Tj7QXpH58/S5Tc3WfuAzsoJk1wXbZuGpRLVpr15rctMeutzZ1VL2pvXmq7Zhyuf8FuxaYy1kzxpvap95uvynavNu+d/sca5ClrbUxyeVJy6MAejw/36ECXuYXkGyG5Sm/AH8yLO/5BTDDooryLu0Qkc8CWI4/hxTg/1E7IorI85yqBPOLKAd8z7Aok639y5ofAeAEFC7gJyh8Unw4pXERURX5fgi+AswvohzwPcOGnGyp6lcBQETuB3C8qu4o3r8awC2pjo6IUldYwsLzpCoT84uo9oWQYXG+bD0YQN/L9nYBmJ7oaIgoE9mc7VFVzC+iGuZ7hsWZbP03gCdE5DYUznd4L4AFqYyKiKrK8w+FSWB+EdUw3zMsznW2vi4i9wB4e7H0UVV9Mp1hWYKhu3Jc3SNNo6Luov1Jbeu0r1cNTaNtV03Uzh2Xjq7+c/7lm+x6Xi51jvf7pGljTe30Q8aVNS4gerdoT8Qv5F1rgbm6mRa+WP77mTbXemHXPdFuav8w23ablaNwQUDPk6pCIeQXUEmG5Se/AH8yjPnlxgyzYvVsqupSAEtTGgsRZcTvmEoG84uodvmeYdlfIIOIMud7Jw8R0WB8zzBOtojI+/MdiIgG43uGcbJFlHMKRa/3B+GJiNxCyDBOtojyTv3/VEhEVFIAGRbsZMu1FtikkXbtPVeHQp1EayH5twfXxB5XXBegc8EAABghSURBVBe1tZpa1M6VqB5dZ9f9GqhxhP1VOLGl0dSWvmLXODt4nF3jbfOebltzrOf1thbbGeQyc5Jdg2zNlmgdSRrx550FV/fRPsfP/+l19nFIqJMH8P98h1qUdoblKb+AbDIs7/kFMMOiCnayRUTJUPj/qZCIqJQQMoyTLSLy/nwHIqLB+J5hnGwRkfefComIBuN7hnGyRZRzIVx9mYiolBAyrC7NFxeRG0Rkk4gs71O7WkReFpFlxdu5aY6BiIagQE9vvFteMMOIAhAzw7KQ9pGtGwFcC+BnA+rXqOq3KnlhVzeOa2Ybda67ZIPtULngTc2m1jjCdgtd+/Dzpnbh8S2mVu9aqMsh6W6eKHbvs+uo3bt6a6TnbnF07bh+Fq79cnWyNNTbzwBnzhhvai927Im0jaj74bMPOrq+khLCp8IM3YgAMizv+QVkk2HMr+jynmGpTrZU9WERmZ7mNoioUooez4MqK8wwohD4n2Gpfo04iEtF5OniIfqJpR4kIvNEZLGILH518+Zqjo8oN/Z/Koxzq4SIXCgiK0SkV0TaktmLqhsyw5hfRNURN8MqUW5+ZTHZ+hGAmQCOBbABwLdLPVBV56tqm6q2TWlqqtb4iPKl+udsLQfwPgAPV/xK2YiUYcwvoiqp7jlbZeVX1bsRVXXj/j+LyE8A3FntMRDRn1X7fAdVXQkA4vmVsUthhhH5pZoZVm5+VX2yJSLNqrqhePe9KMwSiShDvp/v4BNmGJF/fM+wVCdbIvJzAKcBaBKRdgBfAXCaiByLwmT0RQCfiPJaKtJvjSjnnNLxXm/vtMcMx4+w356+pXlMlGFUdI1aV/fR3OnjTM3V4ZK2SjqI9iXcS+t6vfbtdl3F5rF2PbP27XuTHUsF70uDo3sryuvVO365Tz/E/p4kpfCpMPbTmkRkcZ/781V1/v47IrIQwIGO512lqrfH3lpGksqwgfkFJJthec8vwJ8My3t+AUFkWNXzK+1uxA84ytenuU0iiknL+h/VZlUteXKoqs6tbFB+YIYRBSB+hlU9v3gFeaKcU1TeoUNElJUQMiyrSz8QkUd6NN6tEiLy3uJXcicCuEtE7ktiH4gov3zPLx7ZIso51eTPuxt8e3obgNuqtkEiqmnVzLBy84uTLaKcC2GpCyKiUkLIsGAmW6IKKePNdHUe/vNv7VpgXz7z0EivV+/oyBk7ssHUbl32sqldcebhkbZxzswJpramw3azvOBYW6uaRyjS4lqDbf2OTlPb1mnXM6uEq9PG1ZET9blROneidvzc5VgfzbUGW7kqPbROgys3v4BoGcb88kee8qvUNphhVjCTLSJKRwifComISgkhwzjZIso7VfRWcD0eIqJMBZBhnGwR5ZzC/0PwRESlhJBhnGwRkfeH4ImIBuN7hnGyRZRzhU+FfgcVEVEpIWRYMJMtBdDd570cFm/B7X6idu64t2F/oK7uk+bxo0xt5gTb9TNwvTQAzq6lQycOj/R6z23tMrU1W3abms9cyy5s2WO7mbLi+nlHVW6Hoksla8H1o/D+fIfQDcwvIP0MY35lI+/5BTDDXIKZbBFROkI434GIqJQQMoyTLSLy/nwHIqLB+J5hnGwR5ZxCvT/fgYiolBAyjJMtorzTBM+dICKqtgAyjJMtopxT+B9URESlhJBhwUy2du/rxZOv/Lkr5S3NY8xjXA0Qru9x6xwdNHWOLh2XZZvsel5XnD7T1KaMdr21dhuuI58VNCnhcEfXz8vb9/a77+oo8f0X1SfjR9ifrWudM9fvWUN9//uVrAVXSVdRXxrAp8LQDcwvINkMy1N+AcywSiSZXwAzLKpgJltElB7fg4qIaDC+ZxgnW0Q5p1Dvg4qIqJQQMoyTLaK8C+AQPBFRSQFkGCdbRDkXwsmlRESlhJBhdWm+uIjcICKbRGR5n9okEXlARFYV/zsxzTEQ0eD2n1wa55YXzDAi/8XNsCykfWTrRgDXAvhZn9oVAH6rqt8QkSuK97841AuNbqjDcQeOfv1+JVeLdXf3OB7n6KuZPmGEqY0bbls0HI0crmaeikRdl+y06eOHfMyitdtNbW93+V0mvou6xpfrca7OnbNmjDe1FZs7TW3GBNttZZ73qu0Y29drfxaTRtm15cqVpwlUTDcigQwbmF9AshmWp/wq9bg8ZZjP+QUww1xSPbKlqg8D2DKgfD6ABcU/LwBwQZpjIKLB7T+51PdPhllghhH5L26GZSGLc7YOUNUNAKCqG0RkagZjIKIiVaA7RxOoBDDDiDwSQoZ5fYK8iMwDMA8AWltbMx4NUe3K09GqamF+EVWP7xmW6teIJWwUkWYAKP53U6kHqup8VW1T1bampqaqDZAoT3iCfGyRMoz5RVQdIZwgn8Vk6w4AFxf/fDGA2zMYAxH10aMa65ZzzDAiz/ieX6l+jSgiPwdwGoAmEWkH8BUA3wDwSxG5BMBLAC6M+Gr91mqKuubhvzz4vKm1b91taj/5i6NNzbWNrh5bcz2u3tnOY63psB0fsxwdH4vX7zI1VyfHzAm25urcGShqdwsBNy9eZ2qu9cbeepBd+y4K15p5EX+dylLtqy+LyH8AOA9AF4A1AD6qqh1VG0AMyWWYmHxKMsPynl8AMyyqtPMLqO0MKze/Up1sqeoHSvzVmWlul4iiy2AR1wcAXKmq3SLyTQBXIsLlX7LADCPyX5UzrKz88voEeSJKnwLo6u6p3vZU7+9z9zEA76/axomo5lQzw8rNL062iPJOMz3p/WMAfpHVxomoBmSXYZHzi5Mtopwrc12xJhFZ3Of+fFWdv/+OiCwEcKDjeVep6u3Fx1wFoBvATXE3TkS0XxkZVvX84mSLKOfKvCDgZlVtK/2aOnewJ4vIxQDeDeBMVbY3ElH5ysiwqudXMJMtgaKuz+Jc7rXArBNnTDK1WxzdiC7DHNuYNNKuI3bnqoGreQDvnmW363o9V+eOywzHmmaTRrl+fEP/3B9p32lqE0faLqATmkebWlT3rEm2ueycmRNM7f4XtplaNQ4lu7bxRLvd33K7eaJ27UTt0hpKmUe2yiYiZ6NwQumpqhrtH2PgBuYXkH6G1Wp+AeFlWJ7yC6jtDCs3v4KZbBFRSqrfjXgtgBEAHpBCKj+mqp+s5gCIqIZUN8PKyi9OtohyrtrX2VLVw6q2MSKqedXMsHLzi5MtIuISPEQUNN8zjJMtopzL4KKmRESJCSHDONkiIqjnQUVENBjfM6ymJluuDoi5M2wXyOrX7DpdLr2wL/jrla+aWttB403N1eHznsMnm9rADqVSJji6iKJ6bW//HqepY2wH0eETo3UVuSzbtKfs51ZidIN9T3Z0dkd67shhdg32g8ePjPTc3V32SsXbdu8ztVUdXaa2dP32fvf/6g1N5jFJdehEpQr0eh5UeZFkhtVqfgG1kWG1kF8AMyyqmppsEVE5FLzUFRGFy/8M42SLiLw/BE9ENBjfM4yTLaK8C+AQPBFRSQFkGCdbRDmnANR16XIiogCEkGGcbBGR9+c7EBENxvcMC3ay5eq0cS3H5OqW+fu2g6Jtw/HDO6nVdgY1j7Vv4yHj7dpirtdzTcZda5DVRV1symHgL+FhE+06Za41ydSxTVeXyWZHJ0vSNu62HTQntzaaWiVdMM9s6Sz7uVPG2ff0p4+9ZGpbdvXfRsde+94dPtmuSebqXGqoL/93op8ADsHXorQzrFbzCwgvw2o1vwBmWFTBTraIKCnq/cmlRESl+Z9hnGwR5VzhfAe/g4qIqJQQMoyTLaK8U/dXREREQQggwzjZIiLvPxUSEQ3G9wzLbLIlIi8C2AGgB0C3qrZlNRaivPM9qHzEDCPyh+8ZlvWRrdNVdXM5T3R16Cx4apOpNTfaLouzZti1wKJ6euMOu42xE03N1X3jOswZ9dcj6hpkLk2jBnaBRHstV2fM5+9eZWpnHjGlnGHFsnSDfd/r6+x7HPVn+/y28ruPXNtduGJjWa/V3GjXM9vhWLvMVUuKqnrfyeOxoDKsNvIr+pZ9ybBazS+AGRZV1pMtIvKA79eoISIajO8ZZpcOrx4FcL+ILBGReRmOgyj3tDfejQAww4i84Xt+ZXlk62RVXS8iUwE8ICLPqOrDfR9QDLB5ANDa2prFGIlqngZwQUBPDZphzC+i6gghwzI7sqWq64v/3QTgNgCzHY+Zr6ptqto2pamp2kMkyg3t1Vg3GjrDmF9E1eN7fmUy2RKRMSLSuP/PAM4CsDyLsRDlnnKyFRczjMgjMTMsC1l9jXgAgNuk0PEyDMD/qOq9cV7Ata7Yxu17I9XmzrDrg0Xtlrnb0bVxzmF2HbGo7lq1xdSWr99uah87ocXUDhpT3o/P9d65Oo1cXSvfOGeWqT3wfEek7TbU27n9yGG2NnXMcFPbsMOuybV7n+1uWbHZPu7oJtvN1eDYN1U7lmGORqjRw21xx95u+8CUJbesmKKnhydixeRVhuUpvwB/Moz5VZk8ZVgmky1VfR7AMVlsm4gGUP+vUeMbZhiRRwLIMF76gYi8P7mUiGgwvmcYJ1tE5P01aoiIBuN7hnGyRZRzqtU9aVRE/gXA+QB6AWwC8JH9nX1ERHFVM8PKza8sL2pKRJ7o7dVYtwr9h6q+WVWPBXAngC9XvgdElGe+51cwR7YU0q8DZU2H7dqYNNZ2gbxpaqOpubpWHM0dGOaoXXf+kabWXcHhy7c0jzO1R1a/ZmqPrdtmau87crKpLXjarq02sDvmHYc6ujEda6G1b7frb33vf18wNde6YufMtNtweW5rl6nt7LKdMacebH+Ors6d4a4fmsO0xgZTc3V0reqw4+vs9qPr5agpYxN7Le1Nb90ysy3Vvu1qYxB9ib1gDcwvINkMy1N+Af5kGPOrMiFmWLn5Fcxki4hSolrVyRYAiMjXAfwNgG0ATq/qxomotlQ5w8rJL36NSJRzikJQxbkBaBKRxX1u/dYGFJGFIrLccTsfAFT1KlVtBXATgEurv9dEVCviZhgyyC8e2SLKOwW0J/anws2q2lbyJVXnRnyd/wFwF4CvxB0AERGAcjKs6vnFyRZR7lX9EPwsVV1VvPseAM9UbeNEVIOql2Hl5hcnW0R5V/1ztr4hIkeg0Dq9FsAnq7lxIqox1c2wsvIr2MlWvaPzZPwI26FxwkG228HxVPzg8XZT63G0iF4+x67x5WogeaR9p6nNmWbH0jzW/gg+d+oMUxs/wrHIlcPMSaMjbXcgV4fmhJF2m651xa66b7WpnXLIeFMb43ijZk2y63793zrbQeMyqsGecrjesQbZrAm2mynqWnIuXZ5087y223ZalavK3Yh/UbWNeSzJDMt7fgHZZBjzqzIhZli5+RXsZIuIkqFV/hqRiChJIWQYJ1tEeafVPbJFRJSoADKMky2i3FP0eh5URESl+Z9hnGwR5V0GFzUlIkpMABnGyRZRzin8PwRPRFRKCBkWzmRL+nfgTB9vOzSmj2+K9FKupcA+Ndt26bi41hFzdfO8rcV20HRHbCBpHmN/LFGfO6LedrgMfGq01beAURHX6Tp5pl3jbOMuuz7YoeNtp5U43s99PbZbRh2dW67Xe2Wn7eZx6XQ05OzsssWv3fesqY0ebjucdneV9w/9vx9ba2ofnnNIpOdu67TvcVlUy7moKcUhtoMwyQzLU34B/mRY3vMLYIZFFc5ki4hS4/unQiKiwfieYZxsEeVdAOc7EBGVFECGcbJFlHv+BxURUWn+ZxgnW0Q5Vzi51J+rShMRxRFChtmzEatERM4WkWdFZLWIXJHVOIhyr3gIPs6NmGFE3oiZYVnI5MiWiNQDuA7AOwC0A/iDiNyhqn+K+hqutaFcHR+uzkP3mKI9zrWemWMJMuf4XI0xri6dfY5a1O6b4w60a4slqc6x/+ceNtHUnnxlt6m1jrPdNw2OHTv9kHGm5nqPXe/JSRHXUbv450+ZWn2dfcVhjlraDp9sf4au9320Y221cnECFY9vGcb8ii7JDGN+uTHDrKy+RpwNYLWqPg8AInIzgPMBRA4qIkqI+n/1ZQ8xw4h8EUCGZTXZmgZgXZ/77QDemtFYiHJNAe+vUeMhZhiRJ0LIsKwmW65jm+ZAq4jMAzAPAFpbW9MeE1E+aS96u7uyHkVohsww5hdRlQSQYVmdIN8OoG/6tABYP/BBqjpfVdtUta1pypSqDY4ob3iCfGxDZhjzi6h6fM+vrI5s/QHALBGZAeBlABcB+GBGYyHKtwAuCOghZhiRLwLIMNGo7XpJb1jkXADfBVAP4AZV/foQj38VwFoATQA2pz/C1NXCftTCPgC1sR999+EQVY18KEVE7i0+P47Nqnp2zOfUlDgZ1ie/gNr7fQtZLexHLewD8Of9iJVfQFkZVvX8ymyyVS4RWayqbVmPo1K1sB+1sA9AbexHLexDXtTCz6oW9gGojf2ohX0Aamc/SsnsoqZEREREecDJFhEREVGKQpxszc96AAmphf2ohX0AamM/amEf8qIWfla1sA9AbexHLewDUDv74RTcOVtEREREIQnxyBYRERFRMIKabInI2SLyrIisFpErsh5PVCJyg4hsEpHlfWqTROQBEVlV/K9dCdUjItIqIotEZKWIrBCRy4r1YPZDREaKyBMi8lRxH75arM8QkceL+/ALERme9ViHIiL1IvKkiNxZvB/cPuQN8ys7tZBfADMsZMFMtkSkHsB1AM4B8AYAHxCRN2Q7qshuBDDwmh5XAPitqs4C8NvifZ91A/icqh4FYA6ATxXf/5D2oxPAGap6DIBjAZwtInMAfBPANcV92ArgkgzHGNVlAFb2uR/iPuQG8ytztZBfADMsWMFMtgDMBrBaVZ9X1S4ANwM4P+MxRaKqDwPYMqB8PoAFxT8vAHBBVQcVk6puUNWlxT/vQOEfyTQEtB9asLN4t6F4UwBnAPhVse71PgCAiLQAeBeAnxbvCwLbhxxifmWoFvILYIaFLKTJ1jQA6/rcby/WQnWAqm4ACkEAYGrG44lMRKYDOA7A4whsP4qHrpcB2ATgAQBrAHSoanfxISH8Xn0XwBcA9BbvT0Z4+5A3zC9PhJxfADMsVCFNtsRRYytllYnIWAC/BnC5qm7PejxxqWqPqh6LwsLBswEc5XpYdUcVnYi8G8AmVV3St+x4qLf7kFP8GXkg9PwCmGGhymoh6nK0A2jtc78FwPqMxpKEjSLSrKobRKQZhU8pXhORBhSC6iZVvbVYDm4/AEBVO0TkIRTO35ggIsOKn6p8/706GcB7iuvyjQQwDoVPiSHtQx4xvzJWS/kFMMNCE9KRrT8AmFXsWBgO4CIAd2Q8pkrcAeDi4p8vBnB7hmMZUvE79esBrFTV7/T5q2D2Q0SmiMiE4p9HAZiLwrkbiwC8v/gwr/dBVa9U1RZVnY7Cv4EHVfWvEdA+5BTzK0O1kF8AMyxoqhrMDcC5AJ5D4Tvqq7IeT4xx/xzABgD7UPiEewkK31H/FsCq4n8nZT3OIfbhbSgc1n0awLLi7dyQ9gPAmwE8WdyH5QC+XKwfCuAJAKsB3AJgRNZjjbg/pwG4M+R9yNON+ZXpPgSfX8X9YIYFeuMV5ImIiIhSFNLXiERERETB4WSLiIiIKEWcbBERERGliJMtIiIiohRxskVERESUIk62iIiIiFIU0hXkqUpE5GoAO1G4su/Dqrqwgtc6FsCPiq/VA+DrqvqLJMZJRDQQ84t8xMkWlaSqX07gZXYD+BtVXSUiBwFYIiL3qWpHAq9NROTE/CKf8GtEAgCIyFUi8qyILARwRLF2o4i8v/jnF0XkX0XkURFZLCLHi8h9IrJGRD5Z6nVV9TlVXVX883oU1h6bUoVdIqKcYH6R73hkiyAib0FhjarjUPidWApgieOh61T1RBG5BsCNKCwoOhLACgA/jrCd2QCGo7BcCRFRxZhfFAJOtggA3g7gNlXdDQAiUmqB3P31PwIYq6o7AOwQkb0iMmGwQ+si0gzgvwFcrKq9CY6diPKN+UXe49eItF+URTI7i//t7fPn/fdLTtxFZByAuwB8SVUfK3uERERuzC/yGidbBAAPA3iviIwSkUYA5yX1wiIyHMBtAH6mqrck9bpEREXML/Iev0YkqOpSEfkFgGUA1gL43wRf/i8BnAJgsoh8pFj7iKouS3AbRJRTzC8KgahGOfpKREREROXg14hEREREKeLXiJQIEXkTCt06fXWq6luzGA8RUVTML0obv0YkIiIiShG/RiQiIiJKESdbRERERCniZIuIiIgoRZxsEREREaWIky0iIiKiFP1/H6nDJyk2uCQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=[10, 10])\n",
    "ax1 = plt.subplot(221)\n",
    "ax2 = plt.subplot(222)\n",
    "xds_c_read[-1,:,:].plot(ax=ax1)\n",
    "xds_f_read[-1,:,:].plot(ax=ax2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure the fortran arrays are equal\n",
    "np.testing.assert_array_equal(pfb_f_data_read, pfb_f_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The C arrays appear equal until the last Z layer\n",
    "np.testing.assert_array_equal(pfb_c_data_read[:-1, :, :], pfb_c_data[:-1, :, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "\nArrays are not equal\n\nMismatched elements: 95 / 84050 (0.113%)\nMax absolute difference: 3.6163819\nMax relative difference: 15761.57464594\n x: array([[[ 9.800360e+01,  9.800901e+01,  9.800599e+01, ...,\n          8.486828e+01,  9.387037e+01,  9.736460e+01],\n        [ 9.800431e+01,  9.793230e+01,  9.792827e+01, ...,...\n y: array([[[ 9.800360e+01,  9.800901e+01,  9.800599e+01, ...,\n          8.486828e+01,  9.387037e+01,  9.736460e+01],\n        [ 9.800431e+01,  9.793230e+01,  9.792827e+01, ...,...",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-18-eb53f59c46d7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Make sure the C arrays are equal\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtesting\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massert_array_equal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpfb_c_data_read\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpfb_c_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/miniconda3/envs/pfio/lib/python3.7/site-packages/numpy/testing/_private/utils.py\u001b[0m in \u001b[0;36massert_array_equal\u001b[0;34m(x, y, err_msg, verbose)\u001b[0m\n\u001b[1;32m    934\u001b[0m     \u001b[0m__tracebackhide__\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m  \u001b[0;31m# Hide traceback for py.test\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    935\u001b[0m     assert_array_compare(operator.__eq__, x, y, err_msg=err_msg,\n\u001b[0;32m--> 936\u001b[0;31m                          verbose=verbose, header='Arrays are not equal')\n\u001b[0m\u001b[1;32m    937\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    938\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/pfio/lib/python3.7/site-packages/numpy/testing/_private/utils.py\u001b[0m in \u001b[0;36massert_array_compare\u001b[0;34m(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf)\u001b[0m\n\u001b[1;32m    844\u001b[0m                                 \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mheader\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mheader\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    845\u001b[0m                                 names=('x', 'y'), precision=precision)\n\u001b[0;32m--> 846\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mAssertionError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    847\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    848\u001b[0m         \u001b[0;32mimport\u001b[0m \u001b[0mtraceback\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: \nArrays are not equal\n\nMismatched elements: 95 / 84050 (0.113%)\nMax absolute difference: 3.6163819\nMax relative difference: 15761.57464594\n x: array([[[ 9.800360e+01,  9.800901e+01,  9.800599e+01, ...,\n          8.486828e+01,  9.387037e+01,  9.736460e+01],\n        [ 9.800431e+01,  9.793230e+01,  9.792827e+01, ...,...\n y: array([[[ 9.800360e+01,  9.800901e+01,  9.800599e+01, ...,\n          8.486828e+01,  9.387037e+01,  9.736460e+01],\n        [ 9.800431e+01,  9.793230e+01,  9.792827e+01, ...,..."
     ]
    }
   ],
   "source": [
    "# Make sure the C arrays are equal\n",
    "np.testing.assert_array_equal(pfb_c_data_read, pfb_c_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cleanup \n",
    "\n",
    "Because we don't want to run out of memory on subsequent runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del(pfb_f_data)\n",
    "del(pfb_c_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del(xds1)\n",
    "del(xds2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}