LME_MMR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import imp\n",
    "import ens_snapshot_tools as ens\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd\n",
    "# for debugging\n",
    "# imp.reload(ens)\n",
    "\n",
    "# Luke's\n",
    "import utils #utils file with functions\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import xesmf as xe\n",
    "from scipy import stats\n",
    "from scipy import signal\n",
    "import time as tm #we run into conflicts if we call 'time' from the cesm time\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature\n",
    "from cartopy.util import add_cyclic_point\n",
    "%matplotlib inline\n",
    "import cartopy.util as cutil\n",
    "from cartopy.mpl.geoaxes import GeoAxes\n",
    "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
    "from mpl_toolkits.axes_grid1 import AxesGrid\n",
    "from mpl_toolkits.basemap import Basemap, maskoceans\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "import utils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Luke's system-specific file paths here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATAPATH_MMR  =    '/Users/lukeaparsons/Documents/Data/CESM/CESM_LME/MMR_VolcForcing/'\n",
    "\n",
    "SAVEPATH  =    '/Users/lukeaparsons/python/SEOFs/output/'\n",
    "FIGUREPATH  =  '/Users/lukeaparsons/python/SEOFs/figures/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load LME IVI forcing file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_dataset(DATAPATH_MMR + 'IVI2LoadingLatHeight501-2000_L18_c20100518.nc') #load forcing file\n",
    "LAT = ds['lat'][:]\n",
    "lev = ds['lev'][:]\n",
    "\n",
    "ds_MMRVOLC = ds['MMRVOLC'][:,:,:]\n",
    "ds_MMRVOLC_levmean = ds_MMRVOLC.mean(dim='lev')\n",
    "ds_MMRVOLC_levmean_weighted = utils.weightdata(ds,ds_MMRVOLC_levmean) #use function to latitude weight data\n",
    "ds_MMRVOLC_levmean_weighted_areamean = ds_MMRVOLC_levmean_weighted.mean(dim='lat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plot monthly data to check if weighting, averaging messed up data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'MMR IVI')"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAERCAYAAABxZrw0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAApo0lEQVR4nO3deZwcZ33n8c9vTmlGkiVbYxnLloUvEWLABnE6wcFcBgKEwGZNAphA4uyGmGPDErNswATYEFjAQAhZgbGBgI0xNjYIX/IRy4dkS7Jl67JlWYcly9KMztFoNOdv/6jqnp6ePqp7pvqo/r5fL73U013Hr7uqfvXU8zz1lLk7IiKSPE3VDkBEROKhBC8iklBK8CIiCaUELyKSUErwIiIJpQQvIpJQNZfgzexHZrbXzNZN0fK+ZmbrzWyjmX3HzGwqlisiUutqLsED1wAXTcWCzOx1wPnAS4FzgFcCF0zFskVEal3NJXh3vw/Yn/memZ1hZreZ2WozW25mL4q6OGAa0Aa0A63AnikNWESkRtVcgs9jCXCZu78C+DTwb1FmcveHgHuA3eG/2919Y2xRiojUkJZqB1CMmc0AXgf8MqP6vD387E+Bf8ox2y53f6uZnQn8HnBK+P6dZvb68CpBRCTRaj7BE1xlHHT3c7M/cPcbgRsLzPseYIW7HwEws1uB1wBK8CKSeDVfRePuh4GtZvZfACzwsoiz7wAuMLMWM2slaGBVFY2INISaS/Bmdi3wELDIzHaa2UeBvwA+amZrgfXAuyMu7gZgC/AEsBZY6+6/iSFsEZGaYxouWEQkmWquBC8iIlOjphpZ586d6wsXLqx2GCIidWP16tU97t6V67OaSvALFy5k1apV1Q5DRKRumNn2fJ+pikZEJKGU4EVEEkoJXkQkoZTgRUQSSgleRCShlOBFRBJKCV5EJKGU4EUiGB11rl/1LMMjo9UORSQyJXiRCG5YvZPP3PA4P1i+tdqhiESmBC8SwYGjg+P+F6kHSvAiIgmlBC8iklBK8CIl0PMTpJ7EmuDNbLaZ3WBmm8xso5m9Ns71icRl7HnvIvUj7uGCvw3c5u7vM7M2oCPm9YmISCi2BG9ms4DXAx8GcPdBQF0QREQqJM4qmtOBbuBqM3vUzH5oZp3ZE5nZpWa2ysxWdXd3xxiOyOSpCl7qSZwJvgV4OfB9dz8P6AMuz57I3Ze4+2J3X9zVlfOpUyIiUoY4E/xOYKe7rwz/voEg4YvULTW2Sj2JLcG7+/PAs2a2KHzrjcCGuNYnIiLjxd2L5jLgZ2EPmmeAv4x5fSKxUh281JNYE7y7PwYsjnMdIpVgqG5G6o/uZBURSSgleBGRhFKCFymBquClnijBi0Sg7pFSj5TgRSJQ7xmpR0rwIiVQQV7qiRK8iEhCKcGLlEA1NVJPlOBFIlAjq9QjJXgRkYRSghcRSSgleJESqLuk1BMleBGRhFKCFxFJKCV4kRKoN43UEyV4kRKoDl7qiRK8iEhCKcGLiCSUEryISEIpwYuUwDUajdQRJXiRCEzdZ6QOKcGLiCRUS5wLN7NtQC8wAgy7++I41ycSF1f/SKlDsSb40BvcvacC6xERkQyqohGJQHXwUo/iTvAO3GFmq83s0lwTmNmlZrbKzFZ1d3fHHI6ISOOIO8Gf7+4vB94GfMzMXp89gbsvcffF7r64q6sr5nBERBpHrAne3Z8L/98L3AS8Ks71iYjImNgSvJl1mtnM1GvgLcC6uNYnUgnqTCP1JM5eNPOAm8LGqRbg5+5+W4zrE4mNmlilHsWW4N39GeBlcS1fpJJUcJd6pG6SIiIJpQQvIpJQSvAiEagOXuqREryISEIpwYuIJJQSvIhIQinBi0SgscakHinBi4gklBK8SAQaokDqkRK8iEhCKcGLRKA6eKlHSvAiIgmlBC8iklBK8CIlcLW2Sh1RgheJQFXwUo+U4EVEEkoJXkQkoZTgRSJQzbvUIyV4EZGEUoIXiUCNrFKPlOBFRBJKCV5EJKGU4EVKoMZWqSexJ3gzazazR83st3GvSyQuptHGpA5VogT/CWBjBdYjIiIZYk3wZnYK8A7gh3GuR0REJoq7BH8l8BlgNN8EZnapma0ys1Xd3d0xhyMyORprTOpJbAnezP4Y2OvuqwtN5+5L3H2xuy/u6uqKKxwRkYYTZwn+fOBdZrYNuA640Mz+I8b1icROba1ST2JL8O7+WXc/xd0XAhcDd7v7B+Jan4iIjKd+8CIlUB281JOWSqzE3e8F7q3EukTioKoZqUcqwYuIJJQSvIhIQinBi5TANRqN1BEleJEEuHHNTh7Ztr/aYUiNqUgjq4jE639cvxaAbV99R5UjkVqSN8GbWS+5R0c1wN19VmxRidQo07OdpI7kTfDuPrOSgYjUA9XBSz3JWwdvZr81sz83s45KBiRSi1Rul3pUqJH1B8A7ge1m9gsz+xMza6tQXCIiMkl5E7y73+zu7wcWADcClwA7zOxHZvbmSgUoIiLlKdpN0t373f0X7v4e4C3AecBtsUcmIiKTUjTBm9k8M7vMzB4Afg3cAbwi7sBEapEGG5N6Uqib5F8D7wcWEVTRfMbdH6hUYCI1JRxtTPld6kmhG51eB3wVWObueR+5JyIitalQgv9u+P+5lmOsVHdfE0tEIiIyJQol+G8U+MyBC6c4FpGalSriqA5e6kmhO1nfUMlARGqZHvgh9UijSYqUREV4qR9K8CIRpAYZUxWN1BMleJEIVEUj9ajkBG9mi8zsB3EEI1LrVIKXelJoNMmXmtkdZrbOzL4c3tH6K+AuYEPlQhSpvnQvGtXBSx0pNprkz4H3At3AGuAZ4Ex3/1YFYhMRkUko1A++3d2vCV8/aWafBi5395EoCzazacB9QHu4nhvc/QuTCVZERKIrlOCnmdl5jF2dHgFeauFtrRHuZB0ALnT3I2bWCtxvZre6+4pJRy1SYalGVtXBSz0plOCfB76Z5++id7K6uxOcFABaw386PKQu6VmsUo8K3cn6R5NduJk1A6uBM4HvufvKHNNcClwKsGDBgsmuUiRWKqFIPSk0XPCfFprR3W8stvCwvv5cM5sN3GRm57j7uqxplgBLABYvXqzjR2qTCvBShwpV0dwAPBb+g/G7uBOMER+Jux80s3uBi4B1RSYXqVmqg5d6UijBvxf4r8BLgZuBa9396agLNrMuYChM7tOBNwH/MplgRapF/eClHhV66PZN7n4xcAGwBfiGmd1vZhdEXPYLgHvM7HHgEeBOd//tpCMWEZFICpXgU44Bh4DDwAJgWpQFu/vjBA/oFhGRKijUyPoGgmeyvgpYBnzb3VdVKjCRWpJ+qplqaKSOFCrB3wU8DtxPcDfqh8zsQ6kP3f3jMccmUjPUiUbqUaEE/xFUXhEZRwdEfXvuYD8nz55e7TAqptCNTtdUMA6RmjY2VIFSfL16aMs+3v+DFXz74nN597nzqx1ORRSqg7+l0Izu/q6pD0ekNumBH/Vv4+7DADy646ASPPBa4FngWmAlqoYUURWN1JVCCf4k4M0EPWn+HFhKcLPT+koEJlJL9EzW+teI1WyFbnQacffb3P0S4DXA08C9ZnZZxaITEZGyFbzRyczagXcQlOIXAt+hhDFoRESkego1sv4YOAe4Ffhi9iiQIo1E9znVv0ZsRCxUgv8g0AecDXzcxroRGMHzPGbFHJuIyJRrpJN0oX7whR7ILdKQGqmBTuqfkrhIBKkrWKV3qSdK8CIRNGL9bdKkT9INdJZWghcpRQMlh6RpxLuRleBFImjE5CD1TwlepAR6ZJ/UEyV4EWkojXSSVoIXicDUzFr3GnELKsGLlKCRemBI/VOCF4lgbCTC6sYhk9dI21AJXiSCRry8T5wG7AoVW4I3s1PN7B4z22hm683sE3GtS6RSGqmBLqkaaQsWHC54koaBv3f3NWY2E1htZne6+4YY1ykSC1XRSD2KrQTv7rvdfU34uhfYCDTGgxAlgRrv8l7qX0Xq4M1sIXAewbNdReqWCvD1K3WKbqSrsNgTvJnNAH4FfNLdD+f4/FIzW2Vmq7q7u+MOR0QaVAO2scab4M2slSC5/8zdcz7qz92XuPtid1/c1dUVZzgiZWvE5CD1L85eNAZcBWx092/GtR6RSmqky3upf3GW4M8neOzfhWb2WPjv7TGuTyQ2YwX4IMN//fZNfO+ep6sVjkxK45ylY+sm6e73o64HkhCWVUfzvXu2APCxN5xZjXCkDI04npDuZBUpgapopJ4owYtEkO5iV9UoJFv/4Agv+cLtLNuwJ/I8jXSSVoIXiUC9aGrT9v199A4M87XbNxWdthG3oRK8SAm8kYp/dSBVr17KZmmkTagELyJ1q5RSeQMW4JXgRUSSSgleJIL0aJLVDUOkJErwIlL3SjnxNtKY/krwIiVopAa6ejA2QmTxDdOIY/orwYuUoIFyQ11oxK6PpVCCF5G6pxNvbkrwIiVQP/haE731OzWeUCNtQSV4kQjy5fUVz+xj7+FjlQ1G0srpBz/aQCdpJXhpeIPDozz4dE9Z8168ZAVv/879UxyRlCpKyrYG7OuqBC8N76u3buLPf7iSJ3YeKmv+niMDUxyRRFVKG2sjDhiX6AR//+YeBodHqx2G1Linu48AsK+veKJuoKv7ulKsbWTJfVv4+1+ujTRtksT2wI9qW/vsQT5w1Ur+8vyFnNE1g0P9Q3o4g+TUiCW7pIjacPofK3akX4820IZObII/cHQQgC3dfVz9wDZAT9+RIiIc+I10F2Q9iFpFk9kY20hbMLFVNOkzewNdjkl5SumJod1pvP7BkZqoBi1tuODG2YiJTfAiUY1V0TTOgT9Vfu/zt/Gn338g7+d9A8PsOtgf2/qjnpwzJ2ukrZzYBF9PdzAPj4wyNFL9UlCjGrvaKz5tAxX+Ilu363Dez977/Qc5/6t3xx5DKSdnleATpB625av/z1289Io7qh1Gwyu0r9TDflSLNj3fG+vyoz7RyRp00JrENrKO3dNQ+0fmvr7BaofQ0NSLpn6VVUXTQBs6thK8mf3IzPaa2bq41lFw/XVVSSPVNDaMbPEjvx4KDJJDRjrQUAVT4xrgohiXX1Ajjv0s5VJhoJE0Uk6ILcG7+33A/riWLzLVohz3jZQc6knROvjMaWONpLZUvZHVzC41s1Vmtqq7u3vqlhv+rwNSiinlak+7U/1rpJxQ9QTv7kvcfbG7L+7q6pq6BddRI6tU11jpTvtKI1A3yQSI2n1KpKQedNqfalKxpJ3ZTTI15baePhZevpQ1Ow7EGFl1JTbBR7Vu1yE+eNXKmrjdWuLz8Nb9PL03d5/sUgoDk7ki3PDcYT541UoGhkfKXoaMl9pmpWyV1Mlg+eagSvjGNTunOKraEWc3yWuBh4BFZrbTzD4a17oKKbbhP3PD4yzf3MNTe+K9IUOq68/+30O86Zv35fwsynMgpqLg/rlfP8HyzT2s21XeuPOSX1mNrCXcwVyvYrvRyd3fH9eyo7CId6+oO6WkxL0PqOG/sC/+Zj1/eNZcLnzRvJLnHSlaRTP2Onu44CRvjsRW0UQdQKqe7niVeDSFO0GUG2Amk5wb9Xb5qK5+YBsfuWZVSfOkjtuidfCZTenhtI1wwk1ugo94MOmOV0ntKpW6wzHB+aRqytl0YymiOltkz+Fj3LVxT6zrSHCCD/4v9vQWVdFIUymjSU5qPUReT9LE1TUxtdhiJ+dxD/wIJ612T7v3/fuDfPTHpV2xlCqxCb4pYqkstd0baXyKRnBsaITdh6KNQx51X4HJJaqxhNJ4+1rcX7l4QW4sw6e2c7ULd8/u7w/XH18ACU7wYb1q1CK8JMqlP13Na/852jjkpZTgp0Ljpff4v3Oxk3NzRqYbK8HXhjifEZvYBN/clGo4izb9ZH/jd373ft753fsnuZTadPemPXX3QJL7noo+7IWV0shadkTUTkapgnKvkHuPDXHwaP7htFNLLbb4pnE3OmWV4Kt8yo2z9iCxCT61QUeKZPipqhd9Ytchnkhg/+YHt/TwkWtW8c07n6p2KLEpaSyayfSimYJl1Ktix2E+r/zKMs79pzuLTlcsSY5L8Fl18HGWoKNQgp+EqHXwjXnhXNyBviEAtu/rq3Ik5YlSvxmlDn4q6klrpcRYDeUmsWND0a4ciyf4sdfpSWvkiirOE37DJ/iUpJaq3J3rVz2bHophb+8xjg1Fv1U+3Rupvmpo0qKUzsb6wcccTCip+1ohwzH9uKkTb7HFNzdNrKIZW0bheX+1eicLL19asKqoHJVo5E18gi92aZja8OVeQta6pU/s5jM3PM537toMwKu+chcfvvrhyPNHvWEs24NP97Dof9/Kof6hkuabalFO8JWqgx9rF0rmvlZI0c4Ok1VCL5rsRtZi+/Y1D24DYMf+o2UGlyem8H9V0UxCsd8uneATetD1HhsGoLt3IP3eimcmPoflilvW8/4lKya8b2WWbr9799MMDI+yvsrtElEOnrF2mOLT9g8O8w83PF5WLFHbhWrZ6KjzhZvX8Uz3kZLmi+s7p5ZatBdNjtEk0yfciO105X6Fp/f25hy1spSCRbkSm+BTv1mxxJ30EnxTxHrfax7cxkPP7Ms/f4k/T1O4Z1X7Z40Sd3YVzdDIKIeO5r7yeGrPEX6x6tmyYmmpQAn+yefjHTRvS/cRfvzQdi796eqS5ou7AFW0Dj4j06WmTR37RauP0t1oy/sO//lUDwC/WfvcuPebK3DCT2yCT4naup7UBF9uCTylKc/OvftQ/4Qddtx6mVwyOzo4zHu//+CkE1Y5JfhP/3ItL/unOya13lzSCWUkvn3trVfeN+V1xZks6v0lWbKPr9FRZ+HlS6esd1axcHL1oolaZRZx3MLi82ctIPIJZhISm+BTJdbMxsFcZ+Ck14tOtp4v3zgtFy9ZwWXXPpq3f3yUIXgLeXjrflZvP8CXl24ocwmBKAdP9knw5seCE9dU1BuPjDqnf3YpP1+5o2JXi32D8Y03X8pdv5myv3OqRP+vd2+eVDxRw2jKVUVTclfq8rZbvvlbKnDCT2yCT8nceLm2Y0v6oKtURJU12bs00/Nnvb/74DEg/8Ex2frFycZdSmk532iSU1GtcGxohFGHL/12Q8Xae+Js0Cy3x1F2L6yxMWSmIKgImnP0k9x5IBgqoNixb5PcF/NdRTc3B+/HeRNhYhN8rjr4XMlorIom/4/8nn97gK/fvmlqA6wQm2TpIzV/9m9XbATGpkle15YyhG8uYwm++MGT3c7QlOc7lyPzd2oOK4LjLsHH2mhX4nZJ7z9Z00/dvQDRlpNZgk/9/F/53Uag8LEPYyX9cqtS8rWDtaiKZvIyE1uunbI5Qgn+0R0H+d49W6Y8tkpIn8AmWRLOTkrF2i4m27ZRblVASurgGYyS4LOq6SxH7OUegpmjHTZXaFjiUQ9K8Yu/vIzry2wQzieVJ3ce6Gfh5UuLTt+cpwBV6RrRcQX4rK1Z7NhoaZ5kVYrlzjGlFELKldgEn9oUmQdp5utrH97BwsuXMhDeAFQL3SQXXr6Uz9ywdkqXmdo5i5VS8vnINY8AEy8j85XsUyZbHTGWdMuavaT6zeyhpVPJYCpKVmPd+Ma+U7mJYvu+vkhXYqPuDI2O0nNkgP914xNlrSufUsfma8pTgKp0m9e4G52yVl3s2GgNRyobKvMYam3KfQy2hFd0Q6qDL9+4BJ+xZZfc9wwAPUcGwulqoxL++lVT+wDg1M45OOxlVdOkdr7snbBYCbktvd7yftfIo4Hm0RKufzjCdh1rZwjW1TrJ2DOlEtnIqKdLs+Uktyd2HuKCr9+bvumm4DpHvayHUUdRauj5GjJLXc7nb17Hiz9/W9nx5OpFk1LsKrM1Twm+u3eA5w4WH5I6fYLImj/VdTPKPlquxCb4XLcwe8bv2JRVEomzJbua2lrCZDUyGnEwLc/5OrsE39bSDMBAnrFC0uvNSpI7DxzlhtXFT2KZcZcjVWKLUjpK17mH06bWPTA8+d4omftclOrAfFJ3UT6ybeJNatmGRjydtKa6vr/k+yHyVEuVepL7yUPbOTqJ3kFNGQf88OjouH272G+UKmlnV6W88ivLeN1Xiw9J3ZpnX1YJfgpk7kiZJfjUj5tZR5pErU2pRDsSqTSXWS2Rud9nJ/j2cKfNHNemb2CYT/3iMfYdGUiXelLz9R4bYnB4lIuXrODTv1xbNHmmrhAGhkbZ0n2En67YHiH6ifNHOXG3hyer1AGYuvqIOtBVIZn7VXOeS/WUrT19vOJLd7IrR6kwNZ55lO9zbHhk3Hqn8oESxY6TD161kqWP707/3ZSnIXHCg6/LiLH32BD/fGu0zg+ZBbpjQ6Pj1l+0BD/JwkZb6ljIKuyoDn4SitXBp+qmxy6hp2i9NXqiGBqJVkWTmbAzf6/sUsZYgh/74X61Zic3PbqLK5dtzrgsDT5/yRV38OGrH05XiRUrtaRCHRge4X3ff5B//PW6kkqjzSU0sk5rHX+yastx8irX6LhCReGG5+se2cG+vkFueWziDWSpHjhR2gVSXTNTblyzq3ico87PVm4v+p2LrX355h4+9vM16b/z9v2fkODH/90/OMJPH9o2cf0ZE37vni3cvWlv3liODY1w1ud+x2/WPjduqIKBrBNg9q743bs2c+WysRuwWrMKCz1HBvjdE7vHzXPFLeu5ff3zOePIPhZSiu0PUyGxCT4l88fLPNg72sJS2xQ3ssbZ5akcqXrloZHRccfU/Zt70q8z73wcyChlZB4Eg8Oj3L1pD3eEO3E6CWaUxDPrrseqWMaW8eCWfelpjg4GY+Ts7xvklhx3xKbWPTg8yuFwPJ1jQyNs7enj9/7xNpY+vps7N+R/YHFqPdmlo6sf2MrCy5fSe2xsKIJU8kydrMaqaKaiBD/2Ol9p9pa1z/G2by9PH/C5+kV3hvtr6ncrpH9whF8/OpbUN+w+PGGaZ7qPsP65sXGCbl//PJ+7aR3fWlb4ztJCJfhciSrfd8peTvac31r2FP948/oJy3vFl5dx0ZX3hfMUPta6ewcYGnG+euumcY2sQQk+s+A3PrZv3PkUVy4buwEr3YsmnO6vfryKv/3ZmnHzXPPgNv4mz/ANTRndLB/aso+r7t8KZFQj1muCN7OLzOxJM3vazC6Pc13ZUtsv82DKPKhntLcAYwfxwNAId23cE6nRpJCpaJibSqmvPzg8fqf+wFUrcXf6B0d45VeWpd/PLMGNS/Ajo3zkmlXpMUhylXKnt4ZJaGgkbyPrrGmtAGzZ28dt657nE9c9ysevfZTNe3rHnXRS6x4YHk2fjA/1D3HTo7voHxrhYz9fw1//JHhg8dUPbOWXYXfAvoFh3v7t5WztCcavz75SSFX1PBfeqAVjpcJUtVEq9oFJlOB3Heznpyu2jytxppeb9Zt88rpH2bj7cPpEc2RgYhJP/d65PsvWPzTCD5Y/k/47V3XYhd/4T97xnfvD1/fyt2Gpe3fG75JLvvw+NDKas7Fw1vRge2ePKjohwWf9nXmsZtrfN8imcPiKzraWgrFu3xe0WxwdHGbGtLFpB4ZHxn2PYnX7qQb7X60OTprZVWjFSuCpE9Hg8Cjv/8EKvvTbDeOWG2fOKPwLTYKZNQPfA94M7AQeMbNb3H1y955PQuYAUjPDRJM6se/vG+TLSzdy3PRW1n7hLWWv48DRQTrbY/tZS5Y6cI4Ojkw4OPf1DXLdwzvGJcG+gdxVNJmjUQJMC+utM6efNT343gePDvLCuZ3hescnpJnhgfbBq1YyPOrMm9UOwMeve4yNuw9zx6dez9nzZqbLZkcHR+ia2U7vsWG6eweYHSaMlLs27uGLvwl2qTkdbazcum9ciXVf3/i4UwdTrvFaUiNvtocnqlJu+R8aGaW1uYlNzx/mHd+5n462ZnqPDXPuKbPT06SS3YG+YN1/9/M1vGHRicyc1sqh/qGCsaW2xZFjxRN838BweCINElFP79jybl//PNc8sC39t7vzTPfYw1wOFBnHJl8134G+3Pv9nI42oC/9nVOKFVqD+YpN01rw8w9ctTKI7ehQupAAE0vw+44U/s6pq6eHwwbu2dNbxx0PxX6z1HnvcMa2OzY0wgmdwXdMVVvGweKqMzaz1wJXuPtbw78/C+Du/5xvnsWLF/uqVatKXtc7v3v/hLrD/qGR9K3IKTOntXDSrGkAbN9/dNyZ02ysdHJGVydmlh7HZfPeYGjU007oSJfCsqWmmTujveiOl2/elNO7OmnKWP9k9A0M89yhoFR2RlcnWzIO5hM62+gdGB73O8zuaKVrRpB0R7IO/pQzujrZtu8oI6NOR1sz82dPn7CuuTPa6AkPnLNOnDHhO+bT1tLEqXOmj4szM7a25ib29kY/IKa1NnHqnI4J65/d0crcGe0YQQJIHWSZsXa2NXNy+N2ODAyz+1D+0u382dPpaGuesJ5Z01rSB/bxnW3s74v2m5x54oxxf/cPjqRLjmeE+0emzGU1N9mEUuVZ4fKy17ng+I4J45xnrztzTQPDoznHRT9xZjszprWk95fUMrbv62NoxNPbIWVk1HkmvMo668QZOPB0GNupx0/ncP9wwWcJZO/LEBw3mXXtmd91TkcrBzIKeKd3dY7bt8/K+M6p+U6f20lzk7G/b5B94XbLtd4TZ7an98kzT5wx4bjNPC5Suma2c3RgmL7BEaa3NvOS+cdx/X97bd7vW4iZrXb3xbk+i7OoOR/IvI1uJ/Dq7InM7FLgUoAFCxaUtaIzujpzNqad0TWD6a3NNDVBz5FB5s4YKxWcNW8GT+05wqJ5M9l1sJ/5c6azafdhOtpaWHB8x7j6vVPmTOdg/xAvOG5a3hhO7+rk+UPHmD9nesnxnzVvBj29g2zcfZgXnzyLuTPbp7Sx9vj9R5k/ezrNTcbJs6czMDTKpucP8+rTjwfgnk3dtDYbC07o4LTjO8d99zkdbfQeG2JgeJTB4VGazFh00kwWnTQz/ftlTm/PHuLFJ8+itdnYsrePM+fNwN2ZOa2FEYf5s6exZvtBzlswmzs27OGtvz+Puzft5Y0vmsedG/bw5hfPA4dFJ83k2f39nDiznbaWJh7ZdoBXvzCI977N3Zy3YA73PdXN/NnTOXFWO3sPD3Dc9Fb29h5jZNQ5Z/5xHDw6xGkndDDqzuldnRzuH2ZOZyvLN/dw/plzx/3GPUcGecFx0xgaGeWseTN48vleXnTSrHHf7TxgemsLG3cfZsHxHfQODNHZ1sLRwZH01ctZ82bwuyee562/P49lG/fyB2fN5VD/EO7Bb/ncoX5Onj09qB4bGmFGewund3WyfHMPrz+7i1Xb9vPyBXNyPiR6WmsTL5zbSXtL84T657NPmsnwyGhwh+ncTnBYu/NgevkpC47vYPnmHl535gls3H2Yc+bP4oQZbRzuH+LIwDDnnHwc01qbJ6w708mzp9He0syO/Uc5e94Mjg6OpK/MZrS3BPvZccFxsGjeTPb3DXJ8Z9uEZc2d0U5zkzGnMygQnXZ8B32Dw5zQ2Y7jPPD0Ps4/8wR2HzrGycdNZ/1zh5g3axq7DvbzohfMYtFJM3lsx0FefPIsWpqaxg0JDLBwbid3btjDhS86kWmtTfT0DtLZ3kxnewuj7pw4s52eI4MT9uF5s6bx7IGjvOgFM9PvrX32EHNntHHKnA4WnTST5Zt7eM3pJ7C1py+dQ4ZGRjnthA5ymd1zlNNO6GB41Fm9/QCvWng8jrPrQD+nzOlI/35TLc4En6sAOmFvcfclwBIISvDlrOjKi88rZzYRkUSLs5F1J3Bqxt+nAPkHEBcRkSkVZ4J/BDjLzF5oZm3AxcAtMa5PREQyxFZF4+7DZvZ3wO1AM/Ajd5/YsVVERGIRa38+d/8d8Ls41yEiIrkl/k5WEZFGpQQvIpJQSvAiIgmlBC8iklCxDVVQDjPrBkob+Ds+c4GeolNVT63HB4pxKtR6fFD7MdZ6fDC5GE9z965cH9RUgq8lZrYq3/gOtaDW4wPFOBVqPT6o/RhrPT6IL0ZV0YiIJJQSvIhIQinB57ek2gEUUevxgWKcCrUeH9R+jLUeH8QUo+rgRUQSSiV4EZGEUoIXEUmohk3wZjbbzG4ws01mttHMXmtmx5vZnWa2Ofx/Tsb0nw0fHv6kmb21QjF+yszWm9k6M7vWzKZVM0Yz+5GZ7TWzdRnvlRyPmb3CzJ4IP/uOmU3F0wkLxfj1cDs/bmY3mdnsWosx47NPm5mb2dxqxZgvPjO7LIxhvZl9rVrx5YvRzM41sxVm9piZrTKzV1UrRjM71czuCXPLejP7RPh+ZY8Xd2/If8CPgb8KX7cBs4GvAZeH710O/Ev4+sXAWqAdeCGwBWiOOb75wFZgevj39cCHqxkj8Hrg5cC6jPdKjgd4GHgtwVO/bgXeFnOMbwFawtf/Uosxhu+fSjC89nZgbrVizPMbvgFYBrSHf59Ya78hcEdqHcDbgXur+Bu+AHh5+Hom8FQYR0WPl4YswZvZLIId5CoAdx9094PAuwkSP+H/fxK+fjdwnbsPuPtW4GngVcSvBZhuZi1AB8ETsaoWo7vfB+zPerukeMzsBcAsd3/Ig733JxnzxBKju9/h7qlH2q8geLpYTcUY+hbwGcY/2rLiMeaJ778DX3X3gXCavdWKr0CMDswKXx/H2BPkqvEb7nb3NeHrXmAjQaGtosdLQyZ44HSgG7jazB41sx+aWScwz913Q7CBgBPD6XM9QHx+nAG6+y7g/wI7gN3AIXe/o5ZiDJUaz/zwdaXjTPkIQSkIaihGM3sXsMvd12Z9VCsxng38oZmtNLP/NLNX1lh8AJ8Evm5mzxIcO5+thRjNbCHBM9tXUuHjpVETfAvB5d333f08oI/gcimfSA8Qn0ph3dy7CS7XTgY6zewDhWbJ8V41+8Dmi6dqcZrZ54Bh4Gept/LEUtEYzawD+Bzw+Vwf54ml0r9jCzAHeA3wP4Hrw7rgWokPgquMT7n7qcCnCK/QC8QSe4xmNgP4FfBJdz9caNI8sUwqxkZN8DuBne6+Mvz7BoKEvye8JCL8f2/G9JV+gPibgK3u3u3uQ8CNwOtqLEbKiGcnY1UkFYvTzC4B/hj4i/BSt5ZiPIPgRL7WzLaF61tjZifVUIw7gRs98DAwSjBAVq3EB3AJwXEC8EvGqiirEqOZtRIk95+5eyquih4vDZng3f154FkzWxS+9UZgA8FDwS8J37sEuDl8fQtwsZm1m9kLgbMIGj7itAN4jZl1hCWlNxLU49VSjKn1Ro4nvCztNbPXhN/rQxnzxMLMLgL+AXiXux/Nir3qMbr7E+5+orsvdPeFBAf1y8P9tCZiBH4NXAhgZmcTdEzoqaH4IEh8F4SvLwQ2h68rHmO4vKuAje7+zYyPKnu8TEWLcT3+A84FVgGPE+y8c4ATgLsIdoy7gOMzpv8cQcv2k0xhb4AiMX4R2ASsA35K0MJetRiBawnaA4YIktBHy4kHWBx+py3AvxLeUR1jjE8T1G8+Fv7791qLMevzbYS9aKoRY57fsA34j3B9a4ALa+03BP4AWE3QG2Ul8Ioq/oZ/QFCV8njGfvf2Sh8vGqpARCShGrKKRkSkESjBi4gklBK8iEhCKcGLiCSUEryISEIpwUtiWOB+M3tbxnt/Zma3TfE67g7HM8LMTjKz68xsi5ltMLPfmdnZZrbQzPotGNkw9e9D4TzLMkcRFImLuklKopjZOQR3MZ4HNBP0P77I3beUsaxmdx/Jeu8dwJvc/VPhjScPAj92938PPz+XYPTAZ4Hfuvs5OZZ7CXCKu3+l1JhESqEEL4ljwVjlfUBn+P9pwEsIxlO5wt1vDgeA+mk4DcDfufuDZvZHwBcIbqI5191fnLXsnwNL3P1eM7swXN7rc8SwkPwJfg6wPNdnIlOppdoBiMTgiwR3Ww4CvwXudvePWPCgj4fNbBnBGCBvdvdjZnYWwZ2Ri8P5XwWc48GwrdnOB/4mfH0OwZ2T+ZxhZo9l/H2Zuy939wPhLeknuPu+Mr+jSFFK8JI47t5nZr8AjgB/BrzTzD4dfjwNWEAwbsm/hlUqIwTD4aY8nCe5Q3BreW/EULa4+7l5PttLMEqoErzERglekmo0/GfAe939ycwPzewKYA/wMoLOBscyPu4rsNxhM2ty91FgPfC+MuObBvSXOa9IJOpFI0l3O3BZ2CCKmZ0Xvn8csDtM1B8kaJCN4kmCB8YA3A20m9lfpz40s1ea2QU55xybxoCTCAYVE4mNErwk3ZeAVuBxCx7Q/KXw/X8DLjGzFQTVM4VK7ZmWAn8E4EEPhfcAbw67Sa4HrmBsvO4zsrpJfjx8/xXACh97jKBILNSLRqQE4UMafuLub57EMr4N3OLud01dZCITqQQvUgIPHsDwg9SNTmVap+QulaASvIhIQqkELyKSUErwIiIJpQQvIpJQSvAiIgmlBC8iklD/H+shbY1TKfXUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ds_MMRVOLC_levmean_weighted_areamean['time'],ds_MMRVOLC_levmean_weighted_areamean)\n",
    "plt.xlabel('Year (CE)')\n",
    "plt.ylabel('MMR IVI')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "#truncate data to remote dec 500, jan 2001 (at least according to dates in xarray)\n",
    "ds_MMRVOLC_levmean_weighted_areamean_trunc = ds_MMRVOLC_levmean_weighted_areamean[1:18001]\n",
    "time_trunc = ds_MMRVOLC_levmean_weighted_areamean['time'][1:18001]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## now take annual mean of 12 months of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "nyrs = int(ds_MMRVOLC_levmean_weighted_areamean_trunc.shape[0]/12)\n",
    "\n",
    "ds_MMRVOLC_levmean_weighted_areamean_trunc_annual = np.empty([nyrs, 1])\n",
    "time_annual = np.empty([nyrs, 1])\n",
    "\n",
    "index = 0\n",
    "for tt in np.arange(0,nyrs,1):\n",
    "    \n",
    "    OneYear = ds_MMRVOLC_levmean_weighted_areamean_trunc[index:index+12]\n",
    "    ds_MMRVOLC_levmean_weighted_areamean_trunc_annual[tt] = OneYear.mean(dim='time')\n",
    "    \n",
    "    OneYearTime = time_trunc[index:index+12]\n",
    "    time_annual[tt] = int(np.floor(OneYearTime.mean(dim='time')))\n",
    "    \n",
    "    index = index + 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plot annual data to check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'MMR IVI')"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAERCAYAAAB2CKBkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuC0lEQVR4nO3debxcdX3/8dfnLlkJm7kIZQsC2l9LFSRFcAOxVMSFWq3FtsCvVflplaqtbVF/1rUPEatWpBJRqGAVXFiFgBAQSNiTELICSSDLJQm5WW+Su87Mp3+cc+bOnTszd+69c2bO5Lyfj8dNZjlz5jNn+5zvcr7H3B0REUmvlkYHICIijaVEICKSckoEIiIpp0QgIpJySgQiIimnRCAiknJNmQjM7Doz22pmy2s0vyvMbIWZrTKzK83MajFfEZFm0JSJAPgJcG4tZmRmbwTeBLwWOAn4Y+DMWsxbRKQZNGUicPeHgR2Fr5nZ8WZ2j5ktMrP5Zvb71c4OmAJMAiYD7cDLNQ1YRCTBmjIRlHENcKm7nwp8FvhBNR9y98eA3wGbw7/fuvuq2KIUEUmYtkYHUAtmdgDwRuBXBdX7k8P3/hz4aomPveTu7zCzE4D/AxwVvn6fmb01LHWIiOz39otEQFCy2eXuJxe/4e63ALdU+Oz7gMfdfS+Amd0NnA4oEYhIKuwXVUPu3g28aGZ/AWCB11X58Q3AmWbWZmbtBA3FqhoSkdRoykRgZjcCjwGvMbNOM/sw8NfAh83sGWAFcH6Vs/s1sBZYBjwDPOPuv4khbBGRRDINQy0ikm5NWSIQEZHaabrG4pkzZ/qsWbMaHYaISFNZtGjRNnfvKPVe0yWCWbNmsXDhwkaHISLSVMxsfbn3VDUkIpJySgQiIimnRCAiknJKBCIiKadEICKSckoEIiIpp0QgIpJySgQiNbS2ay+Prd3e6DBExqTpLigTSbK3f/shANZd/q4GRyJSPZUIRERSTolARCTllAhERFJOiUBEJOWUCEREUk6JQEQk5ZQIRERSTolARCTllAhERFJOiUBEJOWUCEREUi62RGBmU8zsSTN7xsxWmNlXSkxjZnalma0xs6Vm9vq44hERkdLiHHSuHzjb3feaWTuwwMzudvfHC6Z5J3Bi+PcG4OrwfxERqZPYSgQe2Bs+bQ//vGiy84EbwmkfBw42syPiiklEREaKtY3AzFrNbAmwFbjP3Z8omuRIYGPB887wteL5XGJmC81sYVdXV2zxioikUayJwN2z7n4ycBRwmpmdVDSJlfpYiflc4+6z3X12R0dHDJGKiKRXXXoNufsu4EHg3KK3OoGjC54fBWyqR0wiIhKIs9dQh5kdHD6eCvwJ8GzRZHcAF4W9h04Hdrv75rhiEhGRkeLsNXQEcL2ZtRIknF+6+51m9jEAd58DzAXOA9YAPcDfxhiPiIiUEFsicPelwCklXp9T8NiBT8QVg4iIjE5XFouIpJwSgYhIyikRiIiknBKBiEjKKRGIiKScEoGISMopEYiIpJwSgYhIyikRiIiknBKBiEjKKRGIiKScEoGISMopEYiIpJwSgYhIyikRiIiknBKBiEjKKRGIiKScEoGISMopEYiIpJwSgYhIyikRiIikXGyJwMyONrPfmdkqM1thZp8qMc1ZZrbbzJaEf/8WVzwiIlJaW4zzzgD/5O6LzWwGsMjM7nP3lUXTzXf3d8cYh4iIVBBbicDdN7v74vDxHmAVcGRc3yciIuNTlzYCM5sFnAI8UeLtM8zsGTO728z+sB7xiIjIkDirhgAwswOAm4FPu3t30duLgWPdfa+ZnQfcBpxYYh6XAJcAHHPMMfEGLCKSMrGWCMysnSAJ/Mzdbyl+39273X1v+Hgu0G5mM0tMd427z3b32R0dHXGGLCKSOnH2GjLgWmCVu3+nzDSHh9NhZqeF8WyPKyYRERkpzqqhNwEXAsvMbEn42ueBYwDcfQ7wAeDjZpYBeoEL3N1jjElERIrElgjcfQFgo0xzFXBVXDGIiMjodGWxiEjKKRGIiKScEoGISMopEYiIpJwSgYhIyikRiIiknBKBiEjKKRGIxEDXRUozUSIQiYHygDQTJQIRkZRTIhCJgQoE0kyUCEREUk6JQCQGaiyWZqJEIBIDpQFpJkoEIiIpp0QgEgPVDEkzUSIQEUk5JQKRGLhaCaSJKBGIiKScEoFIDNRGIM1EiUBEJOWUCEREUi62RGBmR5vZ78xslZmtMLNPlZjGzOxKM1tjZkvN7PVxxSNST6oakmbSFuO8M8A/uftiM5sBLDKz+9x9ZcE07wRODP/eAFwd/i8iInUSW4nA3Te7++Lw8R5gFXBk0WTnAzd44HHgYDM7Iq6YROpF3UelmdSljcDMZgGnAE8UvXUksLHgeScjkwVmdomZLTSzhV1dXbHFKVIrqhqSZhJ7IjCzA4CbgU+7e3fx2yU+MmIXcvdr3H22u8/u6OiII0wRkdSKNRGYWTtBEviZu99SYpJO4OiC50cBm+KMSaQeVCCQZhJnryEDrgVWuft3ykx2B3BR2HvodGC3u2+OKyaRNHtpVy+n/fs8NmzvaXQokjBxlgjeBFwInG1mS8K/88zsY2b2sXCaucALwBrgR8DfxxiPSN0k8cY0ty7uZOuefm56akOjQ5GEia37qLsvoHQbQOE0DnwirhhEGiV5aQCCQnoyY5PGKpsIzGwPpbcZIziGHxhbVCIiUjdlE4G7z6hnICL7kwTWDOUlOTZpjLJtBGZ2p5n9lZlNq2dAIhIPq1hRK2lWqbH4R8B7gPVm9gsz+zMzm1SnuESaW4LPunXVsxQrmwjc/XZ3/xBwDHALcDGwwcyuM7Nz6hWgiIjEa9Tuo+7e6+6/cPf3AX9KMFTEPbFHJtLEdNYtzWTURGBmrzSzS83sEeA24F7g1LgDE2lmSWyQtcq9uSXFKnUf/SjwIeA1BFVD/+Luj9QrMBGprXxjcQKTlDRWpQvK3ghcDsxz91yd4hHZLyT5WJvk2KQxKiWC74f/n2wl+p1F9xoQkeagiiEpp1Ii+HaF9xw4u8axiOw3kjjWUCTJsUljVLqy+G31DERkf5LEQ60uKJNy6nKHMhFJDhUIpJgSgUgMkniwVfdRKUeJQCRlEpijpMHGnAjM7DVm9qM4ghHZXyTxyuKojSCJpRVprEqjj77WzO41s+Vm9vXwCuObgfuBlfULUaQJ6WArTWS00Ud/Drwf6AIWE9xW8gR3/24dYhORGCSxtCKNVek6gsnu/pPw8XNm9lngMnfPxh+WSHNL4qG21IWhIlA5EUwxs1MYuiBxL/BaC7cmXVks0pzURiDFKiWCLcB3yjzXlcUiFSTxYKvygJRT6crisyYyYzO7Dng3sNXdTyrx/lnA7cCL4Uu3uPtXJ/KdIiIydpWGof7zSh9091tGmfdPgKuAGypMM9/d3z3KfESaThIbZNVEIOVUqhr6NbAk/IPhJUsnuEdBWe7+sJnNmkBsIk0ryVVDGnROilVKBO8H/hJ4LUEVzo3uvqbG33+GmT0DbAI+6+4rajx/ESmiNCDFKt28/lZ3vwA4E1gLfNvMFpjZmTX67sXAse7+OoJ7H9xWbkIzu8TMFprZwq6urhp9vUh8kniwVfdRKaeaISb6gN1ANzAdmFKLL3b3bnffGz6eC7Sb2cwy017j7rPdfXZHR0ctvl4ktVQzJMUqNRa/jeCexacB84DvufvCWn2xmR0OvOzubmanESSl7bWav0gjJbEeXgUCKadSG8H9wFJgATAZuMjMLoredPd/qDRjM7sROAuYaWadwJeA9vCzc4APAB83swzQC1zgSdx7RMYhyVtyEns0SWNVSgR/xwSqOt39Q6O8fxVB91IRqQMVCKrT3TfI5l19vObwGY0OpW4qXVD2kzrGISJxC+uGklxaSYILfvg4Kzd3s+7ydzU6lLqp1EZwR6UPuvt7ax+OiEhjrdzc3egQ6q5S1dAZwEbgRuAJVLIUqVoSz7rzF5Q1NApJokqJ4HDgHIKeQ38F3EVwUZku+hIZRZIbZJOYpKSxKl1QlnX3e9z9YuB0YA3woJldWrfoRKRm1H1UyqlUIsDMJgPvIigVzAKuZJQxhkQk6WfdiQ5OGqBSY/H1wEnA3cBX3H153aISkZozNfNJGZVKBBcC+4BXA/9QME6JAe7uB8Ycm0jTSvI5d7JLK8nh7qkZn6nSdQTVjEMkIk0iJce0msk5tKZkmelgLxKDJI6WMnQ/goaG0TSSuA7jokQgEoP0HEL2X7kUrUQlApGUiKqGknyNQ5KkaTkpEYjEIMm1CkmOLUnStJyUCERSQt1Hx0aJQEQmKLlHkeRGliyqGhKRCUnk2aQKBGOSyHUYEyUCkZRQ99GxyaVoQSkRiMQgPYeQ/Vea1qESgYhICSkqECgRiMQhyQeRNDWCToSuLBaRCUniwTYtA6iV4u4sf2n3GD8TUzAJpEQgkjYpOsBF/ueJDbz7+wt4+Pmuqj+TpsUUWyIws+vMbKuZlbyPgQWuNLM1ZrbUzF4fVywi9ZbEs8n0lgfg2fCG9Ot39FT9GfUaqo2fAOdWeP+dwInh3yXA1THGIpJ6Q2MNSTVSlAfiSwTu/jCwo8Ik5wM3eOBx4GAzOyKueETqKU0HkWYwntWhxuL6OBLYWPC8M3xtBDO7xMwWmtnCrq7q6/hEGmnrnj6O//xcnt6ws9GhAAUlghQd4IpVUz2WxpJTIxNBqXVSctm7+zXuPtvdZ3d0dMQclsjEOc6ja7aTzTn//ci6Rocj45CmfNnIRNAJHF3w/ChgU4NiEakp9+SdWUajjyYlnqSKzlDVWFwfdwAXhb2HTgd2u/vmBsYjUlNRv/2kHVASFk7iROstTYup7M3rJ8rMbgTOAmaaWSfwJaAdwN3nAHOB84A1QA/wt3HFItII+brPhBxRUnw92bikqS0ltkTg7h8a5X0HPhHX94tIaek5vI1PGkdp1ZXFIjEY3kaQoiNKQo3loD7UuyqeWJJIiUAkBo4PNc4m7ICSpiqPYtVUjw01qqdnOSkRiMQkjWeW+4VwveVStN6UCERi4A4tCasaihJSMqJJvjSVnJQIRGITdR9tcBgyJvnG4oZGUV9KBCIxcJJXNZQvmSQknkaoZl2kcSgOJQKRGLh7wRgqyTqgJKWqqp7Gcg1FUhv546REIBKT/BWqKTqgjNdNT27g0TXbGh0GMJQ00lSlF9sFZSJp5iSvrjnfWJyUgApcdssyANZd/q5Y5j+e35ymkpNKBCJVenZLN6tf3lP19NGZ5ebdfVx281IGMrmYIku2TDbH9+atZk/fYKNDqfI6gkASE2ZclAhEqnTuf87nnO8+XNW0hQeRVZu7uempjTw0hvvlxqFRJYK7lm3mu/Oe54p7nqvvF49TUgcLjJMSgUhMWopOP9PUC6VQVBLaN5BpcCSj6xvM0jeYBVQiSKW+wSxfvmMF3Qkovsr+wEfceqnRjY+e/7++gbQk7cYMFZz0pd+SafSKaoDUNxav3NRN195+Nuzo4SePrsMMvvSeP2x0WNLkfGQeIClHwnqf6Q71wknG76+kMAk0Q7y1kvpEcN6V8wH46vnBwT+TTc/Kl3hZUdVQCk80gaESQbP9/BTlAVUNRdLYQCTxKew+mn+t0VVDYQD1DiMZ/fLH/uVpOhIoEYSG7lPa0DBkP1LcVTGtJxlDF9Y1/vdbiQq7ctK0vpQIQkM7bTJX/u7eQf76x4+zeXdvo0ORKgRtBEW9hhoUS/H3172NoOj7m0W0nPoGs6zt2tvYYGKmRBBqSfhwALc9/RKPrNnO1Q+ubXQoUgV3H1EiSMIZcSNYs2aCMOBP37SEt3/7IXoHsg2OJz5KBKGhqqFkbq0tTdTzQkpr+KrzEQ/qwmjO9reomviRtcEYSAPZ/ffKcCWCUEu+sbjBgZST9PhkGGfkgT9NY9cUyt+gp8l+/oh4myz+sYg1EZjZuWb2nJmtMbPLSrx/lpntNrMl4d+/xRlPRQnfWJt1Z0qz4gN/rsEnlFE8uo6gOlG8YxjBumnFdh2BmbUC/wWcA3QCT5nZHe6+smjS+e7+7rjiqNZQX+dkbqxDY6QnMz4Zzr1UiSAZ6t99dP+4jiCpx4ZaiLNEcBqwxt1fcPcB4Cbg/Bi/b0KSPuKgSgTNZ0TNQkpXXhL2rYkMQ52G+0rEmQiOBDYWPO8MXyt2hpk9Y2Z3m1nJsR3M7BIzW2hmC7u64h3BMak7a4sueGsqjo/Ylhq96oZGH61zY3GCug2N5U5l0WJqoqGSxi3ORFBqkRcvy8XAse7+OuD7wG2lZuTu17j7bHef3dHRUdsoQ9lwrZdrjH1m4y627O6L5burkoirM6VqJdZTtVULO/cNsLtn/xn8sFkv1izOl/vzSViciaATOLrg+VHApsIJ3L3b3feGj+cC7WY2M8aYyhrt8vvz/+sR3nLFA/ULqEhLgq7OlOoUr6lqD4SnfO0+XvfVe2OLp95bUEt4lGm2bTdfNRQ+VyIYn6eAE83sODObBFwA3FE4gZkdbmG50cxOC+PZHmNMZUVdhCut7MEGDkin6wiai+f/KXgtpatu6DqCBgcyRsXx7s/rL7ZeQ+6eMbNPAr8FWoHr3H2FmX0sfH8O8AHg42aWAXqBC7xBpw35A2xCV3azjuCYZiO6jzb4SNKoO5Q1ax17dChKw4CUsQ5DHVb3zC16bU7B46uAq+KMoVqebyNI5spOxgiO+79Nu3r5vYOnTng+Se4+2ijNVzU03P687+nK4lCuQWdL1UrDWUmjzV22mTde/gALVm+b8LyCXkPFL6Zz3SXhV4+r+2jRBWW5BmWC+au7WLR+Z6zfoUQQyuaSXSJoSU4PvP3W8pd2A7BkY212uqSdUeavLK7/FwPJ2LfGcpVwcffRRsV/4bVP8v6rH431O5QIQrlReg01WrMO3JV0P3xoLU+t2wHA1PZWAPoGJz4WRFA1VHwdQTrXXaOGtpiooXibs7F7LJQIQlGJIO6d9aVd47ufgK4sjsc37n6Wv5jzGABTwkTQOxjPcMONPpA06oIyr2GJ4NE125h12V2s2bpnwvMaTX6soRT02FMiCA2G/UfjXNcPPd/Fmy5/gHuWbx7zZ6ONMbsfb4yNNmVSkAh6ajDufM5HXj6W1jVXy95K9658GYD5NWjHGU2ahghRIghF1wjEmfVXbe4G4OkNu8bx6WRfUHbl/av5YHhm3azaw2JXtgbDhOZ8ZGNxo9ddo769lheyTQuT9b7+zND83fnG3FWs2LS7Bt8wZERjcTJ3vZpQIghlwp0/znXdGp7WZyawRSU0D/Cd+57nybCuvVlZldVv1RzQgwJmQldWndXyGp3pk4Me7/sKSm19gzl++PAL+Sq+UsYzcmi9h5hYuamb//jtcw05YVAiCGXyJYL4vqMlf8Y5ni+Jt8SyclN3vttkLud09+0/Y91UKz/U9yjTVbP6srmRJYKk1DHXO4xathFMyTfoDyWC6CDfnylfksuNIxdF0+bbCGK+n8RfXvMYV/1uTU2qJsdKiSAUVQ3FmY1bJ9DoNJ4NeSzOu3I+f3PtEwD84ME1vPbL97Jtb39M3zbkI9c/xRdvWx7791Sj2hJBpUReOBRI0rqPRj+s/uPq165HXqkhraP1UWm9jKedong/HW2/fXZLd74L8nhEs29EO6ASQagejcWtEygRDJ1V1TKi0uYu2wIwYrTVvsEsH5zzGMs6y2/sY73oZt6qrfz08fVjDzIGVuXNiSodEKKhQD590xL+/meLh72XkAJBXfQNZvPbyWglgrGcfJUaRrqaTc7HkQTz1xHk2+cqT3/uf87n3d9fUPX8i0UnEZkGjGmmRBAaaiMYuRJqVUqIqobGVyKoT/dWgLbW0glrxaZunly3gy/eXv4Mvplu8F12WY6yiKN18eiabWzY3jPsvWgdl1oOja4ayjfa1iGMy25eynuuWkDXnv5Rv3cs8QyVCIY+VM3Jx9D+U/13jagainnBRSeKgw3Yh5QIQvleQyXWQa3OwqPG4nGVCKL/67ATRxtkcaN2e5kEUai/6GKsH89/IdbL4xet38mP578wrs+W+xmjLeLo9//Vj5/grd/63bD3Wsdy55MGeXTtdv7xl0ti/Y5lYRXJ7t6BoWqZMtOOZZMuddvLsVSljOVgXnyiUK9EMFChrSMuSgShTLZ+JYLxJPx6DorXHg4gX3zAr+aMpS8zvKHr63etivXy+Pdf/Shfv2vVuD5brg54tPVdKY+3tpRPBJXm+9KuXr43b3WsJb7CWd+y+KXYvgdgUtvQVdpDvYbGXjW0aP0OTv3affkb9ZRqx6mmRBBV2Y1lKPmhqqHwe2Le9aIYG1GqViIIDeavLB75Xq1LBOM5mNdzCOHoRiLlEsFYSgRJVvw7hpJt5c9VOvBUKhBUmu+nbnya7857nme3xH/FbD1Mags2ooFsLn9qNZ4S2PcfWMP2fQP5YUCGhtwqqBqqYp+wfGm8+u2z+J7FcZ+EtalqqHEs30BTvrG4mg3g5e4+vnnPsxUPkuXq3qsRbZRxb4zZnNMWZoJM0U5TTa+G/kz9u77B+Haekd07w9dH+Vyl31+pRFBpvUcnIoVdB2s92mWpM++Xu/v46m9W5rf/WpkcJoI9fRl29wZn8+Uaaitt0lPLDPsxrNdQFftEtFrGUiLI5YKqx2hYmGq3sfGW6qIag8FM/duSYr0fQTNoNSPjXnAdQamqodHn88XblnPvypd58wkzedMJle+2OZ7uYdExOZNz1m/fx7GvmD7meVRjIJMbaiMo2mmiA1mlXg21GLBtPPoGs7S3ju28ZvxVQxUSQYUiQaX5TgkPnP0FB7zBXI7JLa0VY5moL9y6jHmrtnLmazo489W1ux94lAguvu7J/GvlTsYr9eSZOmn4dQOlknU1CTM/RMsYkqsD967Ykn9ebW+ebM7zJ31j0Vqho0HcUl8iyGfhqGqo6P0v3b6cu6sYGyj63N6CS99HTBNONJ4zvegTi9bv5MxvPci6bfvGPI9qDGRy+SJq8QZZrr/2jU9uyD8uLBHUc/z28SSg4oRcbbiVahdaKpUIKiWC6EKpguU3npJj70B2RLffSqKvqHUD5aQSSblsY3GFnzm5rTgRjKzCLVxOz23Zw4sV9o2xlHyyuVw+EUH1JYLxjhzQapWrhuJsP0p9Isgv/ExUNTR8YV//2Hr+8ZfPjDqffBG2wlWB1Vz4Uk5xXOMdxXQ0/dls2d4L0QZeXGVU2Gun8OrOel4Y0zeOEUO9aH8b6mteWTCOUOmpKuSBiolmSnuwKxYmtPEcUD58/VOc/o37S75Xam5t+dJfjauG2kskgnFsD8U91YZuIFXYRjD0+B3/+TBv+48HR8wn+txYlulAJpffr6H6aqXxntGP1hljfCMSVCf1iSB/EUd4cBvvss4PhjVQvkQQbbATaSyOxNXFbCCTyxdri78jirt4g4zO2mD4AXkiG+6P57+QbyCsxngSQfF6qHYo8mzOy/62SlVDlUpIUYmgsI0gO44Lix5dux2o7uw1l/N8ddrjL2wf83dVUrJEUObnVHOBXrQocvl1VN3nAX7x1AZ+vagTGFsi6M/k8vs1VF8iGCzabzbv7mVXz8Conxut++hExigbjRJBPguPPAiM5QwmKkJWKhFEG+x4VmhxPWpc9YhBG8FQj49CUR1p8ZlR1EMEhpcIyv3Oaxe8yGd+sWTYa/mz8fD/r9+1quIgYpHouNs3mOMHD67hsbXVH9DKVQ2Nttpz7mV/W8WqoUqJIEymewvGeCr3Hf/8q2c4q+j6hWJ7+kaekBT/rp7BbP6M+/rHant1d+HJQf77w2140fodzLrsrnwVVmFcxftcdHDM5UsEUaltaLrRdoV/vXlZ/vFYrtodyOaGbdvFJeEzv/U7biqoFo0U7h/ZnHPGNx7gozcsJJPN8b15q8tWH49WIoizN5ESQX5E0JGjj46lh8G0qhLB8P/HovgzsZUIsrn8mEjVlgiGJ4LRz2i/dudKbn16eD/2TM65Z/kWjvvcXDbu6Cn5uVKiax76MlmuuOc5PvSjx6v+bPGZ5FAyCkoY5ZZxfyZXvkRQIRFUOgGIqoYKR9Us9x2/WtTJuu2Vl9GeKgYN7OnP0FZw5l7L6qH2tpHLIVrcP3siOHg+vLoreL1gmuKfXHxxY7ZgHQEsXLeD25dUf01E8cEc4Kl1O3jzNx8YMdDiQCY3bJ0V9uYZzOZYv72Hy25ZRrHogP35W5dxznceCr9jJ3cv38J35z3PFfc8WzK2ocZiZ8Wm3SP2gziHnkh9Imgt6rJV2A9+LBk4qhLYW6FqKF/POYGxhiKVEs5EFJYIHlu7fViR9q5lQaN58c40uSARdPdm+P+3LaNrT3/Jna5Q4Y43kMkxN5x/qbP6e1dsYWnnrhGvR9VYxVVDl928lPdetYA5D60tm1iKw4sSQ38my+9/8R7ee1XpcWP29mfKtn9USgT7KnQkiEoShWfyxcvvfT94hH+/a2X+ealEESWU7t4SJYLieAay+SqTIL6R29Si9TuGbWu/eWYTn/nFknFddBcts+mTgs6K0XwL51X8m1uL2jCKuzB/YM5j/ODBtSO+6/v3r2bxhp0jXi91cve9eavp3NnLwqKqyIFMbtjBd7AgtkojhEbHjZ8/sYEXChquo9Lr1u7Sgzm2FLRXvuvKBbzliuGlvqatGjKzc83sOTNbY2aXlXjfzOzK8P2lZvb6OOMpJd+/OFzJe4oOTsU27+7l/lUvj3g96nVUageMDFUNjaPPe9FuXKl30kQU9hq6Z8UW/jy8Knjhuh38PDyTq9RGcOfSTfzP4xv4ym9WjNpG8NLOoQbvgUyOw2ZMBmDjzqED930rX+bZLd1c8tNFvPeqR1j98h42FTSUR7FGfdUjNz21kaWdu7n87me58Non2NUzwMXXPZlvZL9lcSd/ec3wqqco3GjZlru4a19/pmxpp61CIii1zm57+iUeW7s9v6wKk0Vxsn96wy5+NP/F/POdJeqdo3VRTYmgODEVnxHv7h3k/Vc/xiU/XcjW7j5mf30el974NLc+/RKdOyt3Vii1fPaGSW7a5LAaLPz+win7BsqXwoCSy6mUb9/3fMmqxVLXuSxYEwy/3rVn+AF6IJMbdjJYeDyo1Cb1XJntJlo35doRR7ugbDztYNWK7ToCM2sF/gs4B+gEnjKzO9x9ZcFk7wRODP/eAFwd/l83+aqhcMPtLjgjK7VCLv350yxcv5N7Pv0Wfv/wA/OvR2cslcbxj0oCu3rGPtZ/8TE1rkSws2dwWD33C13BGc0HCnaq4rOqwh4iT60LzsI27+4bdgaTy/mI+vPCg8mu3sH8TUfWF1R7fPSGhcM+c853HwZg3eXvAoJG1u6+DC+XOcsCWLe9h18v6uSh57v40u3LOf6wA/jhQ8PHJ9rTN5hP1IVn5T0DGaZNGr6b7OvPDK8yyOaquoahcJ19855nWbJhF4+FjbQfOu2YEdN07Q0GbPvsr55hzt+cOmJ+2/cOMPOAycNem9Lewu7eytvh0G8bfmDZtrefow+dBsDnblma/03zV2/jiRd3DBuW/OXuvvy0pZQqMUXb/dR8w3iYCAom7drbx0HT2vPPi/eraB0VJ/6SMZQ4Edm2t3yj7fqi6rY9fRlmzhiax/aCz1YqEXz+1mW884+OGPF6dFHcaG0E5e6ZXbh+as3i6ptqZmcAX3b3d4TPPwfg7t8omOaHwIPufmP4/DngLHcv23F/9uzZvnDhwnJvl/XQ8118/c6VI15/Ydu+ERvMCYcdgBEUxcr1SZ55wGQOKdhgt+7pz2+cJxx2QMnP7OoZzO9MJ5aZppxdvYPDzljM4PiOsc0Dhi7RL7Z6614g2EnbW21YQjxu5vQRyyFaRhAc9Ett3LNeMS1fl33EQVPyB/o14XcVipblziqT5Ks6ptNqxovb9o0oMh/fMZ21XWO7zuKVB04mk3W27xt+oJgxuY3DD5oCDC2jKe0tHHHQ1PwyeeWBk5kxJYh/w/aeig350XpfXbQMZkxuY0/RMpwxuY19A5mybUoHT2sfkQiiZTu5rWXEQWNXz8CwA2Hx/GdMaePwA4f/1sgrD5w8LNnOmNLGK8NpS3l5d9+I3wPBdrNj3wA79g3Q2mIcN3M6uZznq1AOnT6JV0yflJ++a29/PoGceNgBbA8/Gz0vjrPU9xVvb8X73mjziJjBCeE+t6W7L3/CcHxHcHHni9v25Zdlqe/tmDE5vw+X2v9f2tVLz0DQfTs6JhVuywdNbefSs0/gI295VVXxjozfFrn77JLvxZgIPgCc6+4fCZ9fCLzB3T9ZMM2dwOXuviB8fj/wr+6+sGhelwCXABxzzDGnrl8/9h4Oi9bv5NoFpUepPGzGFLr29tPeYuEZZmH1kNPearS1ttDeavQNZlm/vSfYgEssu7aWlspDCWRzHDytfVxn9GZGxwGTecNxhzJ3+ZYxtzWMNhZ7fmRHD+6u1TMQ9CqZ1Bb8pt29g5x67CHs7Bkc0R1uMOsce+g05q/eRmuLccJhB5DJ5TAMs5FtHDv2DXDo9ElBZop66+B07enn8IOmsq8/gxFUuQ2GVzs7zuS2Vvb0DXLYjCn535PNOZPaWunuHaS9tYVJbUYm62zb20/OgzP4V79yBi9s28fxHdNZtH4npxxzMBt29HDsK6YzqbUlX2Uwtb2NFgsSUt9glgOnDpUGouW/dU9f/rlR1OsF55hDp7NpVy99g1kOmTaJSW0tZHI5BjJO72Cw3nsHsmze3UfHjMn0Z3J0HDCZwWxwAVMm50xtb82fMS/t3M3rjjqYnoEMWYfpk1rZ1TMYLL8S67h/MJfvjlpsxpQ2JrW1DJ3dGrzhuENZtH4nmaznl+mevgw5dwYzziHT22ltMTbu6OVVHdPZ2t1f8rtLxdLW0kLvYJYp7a1MnzS0b0XrOdqGW1osv38VG8w60ya15kvoU9pawYKqklwOtu/r57AZU5jSHmwbUye1smV3H4cdOCU//97BLCcdeRCdO3pGDIzYYsarZk7nhW37cODgqe0cNLWdTbt6GcjmMDP+6MiDWLmpe1i17p6+DAdPmxQcBzxYlqe/6hU8vX4nfZkgtkwuR4sFx5Vszunc1cuxh06rWD3cYsbOngGyOefQ6ZNoMWNSawt9mSx/+geH82enHDnqsi+lUiKIc4iJUiefxUeiaqbB3a8BroGgRDCeYE499hBOPXZk8bpZlSp6iozXRWfManQI+40LTz+20SGMWZyNxZ3A0QXPjwI2jWMaERGJUZyJ4CngRDM7zswmARcAdxRNcwdwUdh76HRgd6X2ARERqb3YqobcPWNmnwR+C7QC17n7CjP7WPj+HGAucB6wBugB/jaueEREpLRYh6F297kEB/vC1+YUPHbgE3HGICIilaX+ymIRkbRTIhARSTklAhGRlFMiEBFJudiuLI6LmXUBtR08fWJmAtsaHUQFSY8Pkh9j0uMDxVgLSY8PJhbjse5e8sbUTZcIksbMFpa7bDsJkh4fJD/GpMcHirEWkh4fxBejqoZERFJOiUBEJOWUCCbumkYHMIqkxwfJjzHp8YFirIWkxwcxxag2AhGRlFOJQEQk5ZQIRERSTolgFGZ2sJn92syeNbNVZnaGmR1qZveZ2erw/0MKpv+cma0xs+fM7B11iO8zZrbCzJab2Y1mNqXR8ZnZdWa21cyWF7w25pjM7FQzWxa+d6VFt1CLL8Zvhet5qZndamYHNyrGUvEVvPdZM3Mzm9mo+CrFaGaXhnGsMLMrkhajmZ1sZo+b2RIzW2hmpzUqRjM72sx+Fx5bVpjZp8LX67u/BLcl1F+5P+B64CPh40nAwcAVwGXha5cB3wwf/wHwDDAZOA5YC7TGGNuRwIvA1PD5L4H/2+j4gLcCrweWF7w25piAJ4EzCO5kdzfwzphj/FOgLXz8zUbGWCq+8PWjCYZ2Xw/MTOAyfBswD5gcPj8sgTHeG30HwTD4DzZwPR8BvD58PAN4PoyjrvuLSgQVmNmBBBvStQDuPuDuu4DzCRIE4f9/Fj4+H7jJ3fvd/UWC+yycRrzagKlm1gZMI7jDW0Pjc/eHgR1FL48pJjM7AjjQ3R/zYCu/oeAzscTo7ve6e3Qz6ccJ7pjXkBjLLEOA7wL/wvBbuiZmGQIfJ7gPeX84zdYExujAgeHjgxi6K2Ij1vNmd18cPt4DrCI4wavr/qJEUNmrgC7gv83saTP7sZlNB17p4Z3Uwv8PC6c/EthY8PnO8LVYuPtLwH8AG4DNBHd4uzcp8RUZa0xHho+LX6+XvyM4q4KExGhm7wVecvdnit5KRHyhVwNvMbMnzOwhM/vjBMb4aeBbZraRYP/5XBJiNLNZwCnAE9R5f1EiqKyNoFh5tbufAuwjKKaVU6pOLrb+uWG94fkERcTfA6ab2d9U+kiJ1xrdf7hcTA2L1cy+AGSAn0UvlYmlbjGa2TTgC8C/lXq7TByNWIZtwCHA6cA/A78M66qTFOPHgc+4+9HAZwhL/BViiT1GMzsAuBn4tLt3V5q0TCwTilGJoLJOoNPdnwif/5ogMbwcFsUI/99aMP3RBZ8/iqFiZxz+BHjR3bvcfRC4BXhjguIrNNaYOhmqmil8PVZmdjHwbuCvwyJ2UmI8niDhP2Nm68LvWmxmhyckvkgncIsHngRyBAOlJSnGiwn2FYBfMVQ92pAYzaydIAn8zN2juOq6vygRVODuW4CNZvaa8KW3AyuBOwg2JsL/bw8f3wFcYGaTzew44ESCBpy4bABON7Np4VnX2wnqGJMSX6ExxRQWh/eY2enhb7uo4DOxMLNzgX8F3uvuPUWxNzRGd1/m7oe5+yx3n0Ww478+3EYbHl+B24CzAczs1QQdLLYlLMZNwJnh47OB1eHjuscYzu9aYJW7f6fgrfruL7Vo+d6f/4CTgYXAUoKN/BDgFcD9BBvQ/cChBdN/gaAl/zlq2PuhQnxfAZ4FlgM/JehN0ND4gBsJ2iwGCQ5YHx5PTMDs8HetBa4ivBI+xhjXENS/Lgn/5jQqxlLxFb2/jrDXUMKW4STgf8LvXAycncAY3wwsIuh98wRwagPX85sJqnCWFmx359V7f9EQEyIiKaeqIRGRlFMiEBFJOSUCEZGUUyIQEUk5JQIRkZRTIpDUscACM3tnwWsfNLN7avwdD4TjVWFmh5vZTWa21sxWmtlcM3u1mc0ys14LRsKM/i4KPzOvcNRJkbio+6ikkpmdRHBV6SlAK0H/7XPdfe045tXq7tmi194F/Im7fya8wOdR4Hp3nxO+fzLBaJMbgTvd/aQS870YOMrd/32sMYmMhRKBpJYFY+XvA6aH/x8L/BHBeDlfdvfbw4HAfhpOA/BJd3/UzM4CvkRwsdLJ7v4HRfP+OXCNuz9oZmeH83triRhmUT4RHALML/WeSC21NToAkQb6CsHVrwPAncAD7v53FtyQ5kkzm0cwxss57t5nZicSXKk6O/z8acBJHgwHXOxNwP8LH59EcCVrOceb2ZKC55e6+3x33xkOJfAKd98+zt8oMiolAkktd99nZr8A9gIfBN5jZp8N354CHEMwLs1VYVVOlmCY5ciTZZIABEMC7KkylLXufnKZ97YSjCyrRCCxUSKQtMuFfwa8392fK3zTzL4MvAy8jqBzRV/B2/sqzDdjZi3ungNWAB8YZ3xTgN5xflakKuo1JBL4LXBp2LCLmZ0Svn4QsDk8oF9I0LBcjecIbmwE8AAw2cw+Gr1pZn9sZmeW/OTQNAYcTjDAnEhslAhEAl8D2oGlFtzo/Gvh6z8ALjazxwmqhSqVAgrdBZwF4EGPjPcB54TdR1cAX2ZovPjji7qP/kP4+qnA4z50+0yRWKjXkEgMwpuJ3ODu50xgHt8D7nD3+2sXmchIKhGIxMCDG4X8KLqgbJyWKwlIPahEICKScioRiIiknBKBiEjKKRGIiKScEoGISMopEYiIpNz/Anfaa/oZgBekAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time_annual,ds_MMRVOLC_levmean_weighted_areamean_trunc_annual)\n",
    "plt.xlabel('Year (CE)')\n",
    "plt.ylabel('MMR IVI')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save annualized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savez(SAVEPATH + 'IVI2Loading_Annual',\n",
    "         time_annual = time_annual,\n",
    "         MMRVOLC_ANN = ds_MMRVOLC_levmean_weighted_areamean_trunc_annual)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}