Hovmoller_temp_depth.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "netcdf_index loaded.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-18.04/lib/python3.6/site-packages/cmocean/tools.py:76: MatplotlibDeprecationWarning: The is_string_like function was deprecated in version 2.1.\n",
      "  if not mpl.cbook.is_string_like(rgbin[0]):\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import cosima_cookbook as cc\n",
    "import cosima_cookbook.netcdf_index\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import IPython.display\n",
    "import cmocean as cm\n",
    "import pandas as pd\n",
    "\n",
    "from dask.distributed import Client\n",
    "from collections import OrderedDict\n",
    "\n",
    "import sys, os\n",
    "#sys.path.append(os.path.join(os.getcwd(), '..'))  # so we can import ../exptdata\n",
    "#import exptdata\n",
    "#print('Available exptdata keys: ', [k for k in exptdata.exptdict.keys()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cc.netcdf_index.database_url = 'sqlite:////g/data/v45/ahg157/cosima-cookbook.db'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "basedir = '/g/data3/hh5/tmp/cosima/'\n",
    "exptdict = OrderedDict([\n",
    "    ('1deg',   {'model':'access-om2',     'expt':'1deg_jra55v13_iaf_spinup1_A',\n",
    "                'desc': 'ACCESS-OM2','n_files':-12,\n",
    "                'time_units':'days since 1718-01-01','offset':-87658}),\n",
    "    ('025deg', {'model':'access-om2-025', 'expt':'025deg_jra55v13_iaf_gmredi',\n",
    "                    'desc': 'ACCESS-OM2-025','n_files':-30,\n",
    "                    'time_units':'days since 1718-01-01','offset':-87658}),\n",
    "    ('01deg',  {'model':'access-om2-01',  'expt':'01deg_jra55v13_iaf',\n",
    "                'desc': 'ACCESS-OM2-01','n_files':None,\n",
    "                    'time_units':None,'offset':None})\n",
    "])\n",
    "for k in exptdict.keys():\n",
    "    if not('exptdir' in exptdict[k]):\n",
    "        exptdict[k]['exptdir'] = os.path.join(os.path.join(\n",
    "            basedir, \n",
    "            exptdict[k]['model']),\n",
    "exptdict[k]['expt' ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Client</h3>\n",
       "<ul>\n",
       "  <li><b>Scheduler: </b>tcp://10.0.64.25:8786\n",
       "  <li><b>Dashboard: </b><a href='http://10.0.64.25:44466/status' target='_blank'>http://10.0.64.25:44466/status</a>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Cluster</h3>\n",
       "<ul>\n",
       "  <li><b>Workers: </b>6</li>\n",
       "  <li><b>Cores: </b>6</li>\n",
       "  <li><b>Memory: </b>18.00 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: scheduler='tcp://10.0.64.25:8786' processes=6 cores=6>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "client = Client('tcp://10.0.64.25:8786', local_dir='/local/g40/ahg157')\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0.1° Hovmoller\n",
    "\n",
    "This is the problem case, where we keep running out of memory! My suspicions:\n",
    "\n",
    "* inconsistent chunking between temp, temp_WOA and area_t\n",
    "* save temp_WOA to disk so we don't have to keep it in memory\n",
    "* combining both dimensions in sums?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define common chunking for all variables\n",
    "# these chunks are ~3MB -- maybe they could be bigger to reduce\n",
    "# the number of tasks that the scheduler has to deal with\n",
    "chunks={'st_ocean': 7, 'xt_ocean': 400, 'yt_ocean': 300}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using database sqlite:////g/data/v45/ahg157/cosima-cookbook.db\n"
     ]
    }
   ],
   "source": [
    "ekey='01deg'\n",
    "expt = exptdict[ekey]['expt']\n",
    "n_files = exptdict[ekey]['n_files']\n",
    "time_units = exptdict[ekey]['time_units']\n",
    "offset = exptdict[ekey]['offset']\n",
    "\n",
    "temp = cc.netcdf_index.get_nc_variable(expt, 'ocean.nc', 'temp', n=n_files,\n",
    "                        time_units=time_units, offset=offset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Averaging WOA data\n",
    "The strategy here is to load the WOA data with its native chunking, rename the dimensions so that it plays nicely with other xarray stuff (using `.values` means carrying the entire ~5GB in memory -- not sure how the scheduler handled distributing this amongst the workers, maybe they all had a copy or data was distributed along with tasks, which is memory-inefficient somewhere)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using database sqlite:////g/data/v45/ahg157/cosima-cookbook.db\n"
     ]
    }
   ],
   "source": [
    "temp_WOA13 = cc.netcdf_index.get_nc_variable('woa13/01', 'woa13_ts_??_mom01.nc', 'temp').mean('time')\n",
    "\n",
    "# hacky renaming of variables so we can use as a dask array -- there's probably an elegant way to do this\n",
    "temp_WOA13 = temp_WOA13.rename({'GRID_Y_T': 'yt_ocean', 'GRID_X_T': 'xt_ocean', 'ZT': 'st_ocean'})\n",
    "temp_WOA13['st_ocean'] = temp.st_ocean\n",
    "temp_WOA13['xt_ocean'] = temp.xt_ocean\n",
    "temp_WOA13['yt_ocean'] = temp.yt_ocean\n",
    "\n",
    "# force calculation and save to disk -- we could probably skip the .load() and write directly\n",
    "temp_WOA13.load()\n",
    "temp_WOA13.to_netcdf('/local/g40/ahg157/woa13.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can re-open the dataset from disk with the chunking we want. xarray saves the dataset with contiguous chunking, but ideally we could choose the chunking to save it with. Unfortunately, chunking in the x/y dimensions here means we're taking strides off disk which aren't hugely efficient and gives us extra \"transpose\" steps in the task graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rechunk to match ocean.nc data\n",
    "temp_WOA13 = xr.open_dataset('/local/g40/ahg157/woa13.nc', chunks=chunks).temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One last thing, we'll save half our memory usage on this dataset by casting this down to 32-bit floats to match the `temp` data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_WOA13 = temp_WOA13.astype(np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grid metrics\n",
    "Load up the grid areas with the same chunking as we are using elsewhere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using database sqlite:////g/data/v45/ahg157/cosima-cookbook.db\n"
     ]
    }
   ],
   "source": [
    "area_t = cc.netcdf_index.get_nc_variable(expt,'ocean_grid.nc','area_t',n=1,chunks={'xt_ocean': 400, 'yt_ocean': 300})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# I can't really see a better way around this, but it works fine\n",
    "mask = temp.isel(time=0).copy()\n",
    "mask = mask/mask                 ## This seems pretty dodgy to me, but it works!!\n",
    "area = mask*area_t\n",
    "area_sum = area.sum(['yt_ocean', 'xt_ocean']).load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I just saved this to disk because recalculating it every time the kernel crashed or something was killed was killing me!\n",
    "area_sum.to_netcdf('/local/g40/ahg157/area_sum.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The calculation itself!\n",
    "Now the final calculation is using variables which all have the same chunking, and so can be pretty efficiently read off disk by the workers when they need, instead of keeping things in memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_anom = temp - temp_WOA13 - 273.15\n",
    "var = area_t*temp_anom\n",
    "temp_hov = var.sum(['yt_ocean','xt_ocean'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (time: 270, st_ocean: 75)>\n",
       "array([[-6.444772e+13, -6.359090e+13, -6.192200e+13, ..., -9.067360e+11,\n",
       "        -4.290690e+11, -4.316344e+11],\n",
       "       [-1.884902e+13, -1.779290e+13, -1.595039e+13, ..., -9.141511e+11,\n",
       "        -4.316221e+11, -4.323859e+11],\n",
       "       [-4.161272e+13, -4.065789e+13, -3.885833e+13, ..., -9.208673e+11,\n",
       "        -4.346274e+11, -4.339891e+11],\n",
       "       ...,\n",
       "       [-2.193209e+13, -2.107515e+13, -1.919958e+13, ..., -2.350495e+12,\n",
       "        -1.390042e+12, -9.858980e+11],\n",
       "       [-6.366369e+13, -6.217713e+13, -6.014291e+13, ..., -2.355866e+12,\n",
       "        -1.393905e+12, -9.881095e+11],\n",
       "       [-5.797953e+13, -5.663521e+13, -5.498384e+13, ..., -2.363123e+12,\n",
       "        -1.397597e+12, -9.901801e+11]], dtype=float32)\n",
       "Coordinates:\n",
       "  * st_ocean  (st_ocean) float64 0.5413 1.681 2.94 4.332 5.869 7.569 9.447 ...\n",
       "  * time      (time) datetime64[ns] 1985-01-14T12:00:00 1985-02-13 ..."
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# you could wrap this in a progress bar -- took a few hours I think\n",
    "temp_hov.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# save the precious result!\n",
    "temp_hov.to_netcdf('/local/g40/ahg157/temp_hov.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#temp_hov = xr.open_dataset('/local/g40/ahg157/temp_hov.nc')\n",
    "temp_hov_01deg = temp_hov/area_sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'°C')"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA34AAAFuCAYAAADeYMqmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xu8JWdd5/vPL92de1+SdJPsdHdIRoISERC2CQ6jAgEFdIwXIhGFwARznIERhxnk4jmgjMzg6AjOEfG0hDPBwYlMRMiM0RgCqBwh0uEiJjGTGC690xvSSacvJIS+/c4fq6pTu3bVuu+9V6/+vF+v/dprVT1PPU/Vqt1d3/XUJTITSZIkSdL0OmGlOyBJkiRJWloGP0mSJEmacgY/SZIkSZpyBj9JkiRJmnIGP0mSJEmacgY/SZIkSZpyBj9JmiIR8R8j4hf7KPehiHjBcvRJkiStPIOfJE2JiNgEvBz4f/oo/g7g7X0s84KIOBIRv9swLyLiFyLi7yPi4YiYi4j/ERHfVSlzcUTcGBF7ImJ3RPxtRLyymPfsYtnfqP18bzH/OyPiLyLioaL+bRHxosqy3xwRXyrqzEXEH/VYly0R8YGIeLDo799GxI/UymREfD0iVlemrY6I+yMiK9N+MyLujoj9EfEPEfHyHm2fFBHvi4h9EfG1iHhdZd6JEXF9RHy5aP/Z3ZYlSdIwDH6SND1eAdyYmd/sVTAz/xZYFxGzPYq+HHgIuCIiTqrN+23gtcAvAGcCTwQ+DPwwQBHgPgb8JfAE4CzgXwIvrCxjZ2aeXvv5VDHvfwI3A2cDjyva2Vcs+0rgZcDzMvN0YBa4pW0lIuJM4JPAAeA7gY3AO4E/jIgX14rvqfXxRcU2qHoY+OfAeuBK4Lcj4p+2tQ/8CnAh8HjgOcAv1UZcPwn8LPC1LsuQJGlokZm9S0mSJl5EfAx4X2b+t8q0y4BfBf4JsAt4dWb+eTHv94G5zPzVLsv8R+A36QSXV2fm9cX0C4F/AL63CJFNdT8JfCEzX90y/9nAf8vMLQ3zNhb9PSMz9zTM/x3gUGb2PK21KP/vgR8HnpKZRyrT3wD8K+D8zMxiVO//Ap6amZcXZa4HPgf8WmZGy/JvAP4yM/9zy/z7gFdm5l9U+nNhZl5RKzcH/GxmfqKf9ZIkqV+O+EnS9Pgu4K7yTURcDLwfeD2wAfh+4MuV8ncCT21bWER8H7AFuA74IJ3Rv9KldEJjW+g7Ffhe4Poh1gPgQeAe4L9FxI9FxNm1+Z8GXh4Rr4+I2YhY1WN5zwf+uBr6Ch8EzqMzWln6MPD9EbEhIjYA3wd8pG3BEXEK8D3A7S3zzwDOBb5QmfwFOiOPkiQtC4OfJE2PDcD+yvur6IwA3pyZRzLzvsz8h8r8/UWdNlcCf5aZDwF/CLwwIh5XzDsLmO9S9ww6/8d0KwNwbnH9XvXntOycjvIcOkH1PwPzEfFXxUgjxajmvwZ+iM6ppPdHxBu7tLOxpS/zlfmlR+mcZvoS4ArghmJam9+jE+Ruapl/evF7b2XaXmBtl2VKkjRWBj9Jmh4PsTBMbAX+sUv5tXSuZ1ukGMW6HPgAQHHd3VeBlxZFHgRmevTlSI8y0LnGb0Pt5+GizbnMfE1mfhuda+MepjOCSTH/A5n5PDrh9eeBt0XED0XE91VuFFOOwj3Q0peZyvyq99MZ4Xx5tc26iPgN4MnATxVhlYj4vUr7bwa+URRfV6m6joUhXZKkJWXwk6Tp8XcsPGVxB/BtXco/iYWnH1b9OJ1w8rvFXSi/BmzmsdM9bwG2tN0cJjMfAT4F/GT/3W+XmTuAd9MJWfV5BzPzf9BZ/ydn5l9XbhRTnk75UeAnI6L+/95P0dlO/7s2/a/phMKz6dx4ZZGI+FU6N4H5wczcV+nPz1fa/w/FiOk8C0+rfSotp4ZKkrQUDH6SND1uBH6g8v4a4JURcWlEnBARmyPiOyrzfwD4s5ZlXQm8j851g08rfp4FPC0ivisz7wZ+F/jvxWMZToyIkyPiisopl78EvKK4Du8sgIh4akRc12tFIuKMiPjViHhC0feNwL+gc20fEfGKiPjhiFhbzH8hnWvmbm1Z5DvpBNlrIuKcoq8/Dfwy8PpytK5UvP/nwI/W5xXtv4nO6OfzM/PBXutDZ9Tw/yzW6zuAnwP+a2V5J0XEycXbcls23khGkqRhGPwkaXq8H3hRcZpm+ciGV9IJPXvpXAv3eICI+B7g4aabs0TEZjo3b3lXZn6t8nMb8Od0QiF0Hq/wO3RG4vbQOa30x+lcH0dm/g3w3OLn3ojYDWyjE1BL58bi5/j9JJ3HLpxPZ6RuH/D3wLfoPLKCYtqb6Zx+ugf4T8C/zMzG0bkinP0z4GTgDjqnqr4OeFlmNj7/LzNvz8y2Ubn/QOemMHfXTuts89Zi+3yFzufwG+XdVQt3Ad+kM6p6U/H68V2WJ0nSQHycgyRNkYj4D8D9mfmuHuX+GLgmM2/sVk6SJE0Hg58kSZIkTbllP9UzIl4QEXdFxD09br0tSZIkSRqDZR3xKx6w+7/pPEh3DvgM8NOZeceydUKSJEmSjjPLPeJ3MXBPZt6bmQeA64DLlrkPkiRJknRcWe7gt5nO85JKc8U0SZIkSdISWb3M7TU9k2jBuaYRcTVwNcBpp532jO/4ju9oqCJJWkqHDx06+vrg4cMcPHzk6Ps1q04oph/h4YOHAThtzSpOPXHNgmWU9crypUcOdJZ96omrWbNqFUcOHSQPHSRWr+GE1Ws4ePjw0ekAJ6xewyMHDhZ11nD40UfIQweI1ScuqLOw7SMcPHyEA0eyZ9/WrFp1dFq1b6vyCFn0YZC+lWL1mqNlyvLVOgcPH+naTmdbt/etqledNatO4NQT1yzod9s2KLfdoNugXqe6/k19K9up961ap209V+URjjz68IJ55f4ALOjz4Tjh6P5bb6ec17QfVPedXtsNOLoNuq1PfbtVl1Xd1tXtDQxU53CcsGi7Vdezn8+0W51qvfJvt76t29op1f/mHzlwsGvfAE4+6SSkbm677bYHMnPTSvdjki138JsDtlbebwF2Vgtk5jY6z3lidnY2t2/fvny9kyQBsO/BB/jmrnlO2TQDwGe/8jXu2/cIm9edysz60wCY3/sw9+3rBJ3N607l6Y8/h2/umgfglE0zzO/Zz/zeh4+WL23fsQuA2a2b2HDwG3xr12P/DWy46BnM79nPzIa1jct6+uPPYc8dtx2tc9Kmc4/WqSr7Nrf/US7ZfCZPf/w5AEfXqdq3mQ1rO3X27G/t20mbzj1ar963z37lawCNfSvrlW13q1MtX67PoH0rzWxYe3Qd79v3yNE6VXvWnM783k54KrdDU53q53PSpnPZs+b0xm3QVKfb+pTt1Ped6mdebaecXm7LfXds58ADnXZO3Hju0e0HLGi/up5N7Xz2K187uo+W22D7jl0L9p22/bpcH1j4N1J+rqVh9uu2fafafj+fT/WzHqZOt759c9d847Zua6e+PvXPoO3fA4Bv/7YLkLqJiK+sdB8m3XIHv88AF0bEBcB9wBXAS5e5D5KkPpVBaWb9ady375GjQW8Q23fsWhAYN687dVGZMvh8c9c8rDm9c5BcHCSW6gGyrjy4b1M92FzQt0o43bzu1EXrWD0QZc3pfPYrX+Ppj3+sbzPrT+vZdn19murU26nPb+pbWa88oO7Vj6qZDWsXhPd6W/U2ylBVrg9rTgdgTy10j9K3tv2g6bOrOvDAY/07pbKdT9k007V/pXIfrfZxy9qTe9dp+GKjm07AOX1BYGxS7gtHQ/Oa9jr1cAWL97c9e/Yv6Gc1rLX1rV6nqW971pzOzEXPWNBOL/X1mS/a6QTPhV8aHN0P9j7Mt/dcsqReljX4ZeahiHgNcBOwCnhfZt6+nH2QJA2mPPgCmNv/KHP7H+15UNxU55KWsuWBXjVczO99uHG0sFtb5QhfvW/37XuEmT372VApP7Nphvm9D/fsW72PnPtEAP7nF+9ldmv3M4oWHLz3uT7VbVAG7aagXC9f3wbV9anWbzrgB7j1vt1saWmr7fOp962tn9W+3Xrfbi7ZfGaXLfBYOzMb1vLZr3yN7Tt2Mbt1U2vwObBrnhM3zRytN79nP/OPBDPrT+OUri0tdOt9uxv36/mWfedonR6fUVmv/HvYvmMXs5V5M+tP67Sx6Vy+tWsnBx7YyYkbF34+9dDcth99a9dOZi56xoLt1k11W8/vfXhR3/rxP79474Ivdno5OqLco1z9ywZJo1v25/hl5o2Z+cTM/LbMfPtyty9JGtzMhrVsXncqO/Yf5Nb7vsXc/kf7qlPq1Nvds863du1kZsPaoyGhaZSo7WBwbv+jXft2Sm0Esd636mmrbcq+ze1/lO07di04VbKXW+/bvajO/J79i9an3AYAn9q5f8GIXH1UqixfHfWprw90DrZP2TSz4FTIb+6aZ2b9aUe3W72tUrV/1c+nU2/hZ1q206apb8CiU3Wrjm7rPkaTyj7ct++RTtguTzEtT09uWcbc/kfZsf8gn9q5cH7Zz3KdqutWrdO07zTtp/fte4RP7dzPh+/Y0XMEtPxM5/c+vGjfqStPea3uB3P7H21tp/6Zln1rq9O0j35z1/zRL0/qf6tN+0G9b/X1afpsqvubpNEte/CTJE2++kFYGRK2rl3D/F4WBaz79j3C/J79iw72ygPunfuOsGP/wUUHcNUwUlUeQC8KFkUbTQeu5WhNtW9b1p7M3P5HGw9+631b0E7Lges3d80fDcD1kFCGuHrfqqcpNoXSpvWp1mkKZE3hqh6063oFp537jizqX68Q1/SZNoXZpr6VI4bVcFUPpqVqv8rp9TbK8FOWb/rioOlznVl/GlvWnszOfUcat0EZOpu2RbVOuR6tfwu1bdAWzKujff3uO0112j6fpvWo7teL+lbsN91G38o61VHt+T372VA5DbS+Pm37dbd/DySNxuAnSVqkbeRmy9qTmVkP83ub61TDRXmQt2XtyZy77oRFB9Vl+abQUx487th/8LERsso1QfWAVT2ont/bqVdt69b7dvPZr3xt0TqVYbHpgL8pwJS2ru3cobDfA9fq+nQbVauGl/Jgd+e+I4vKlwfV1e1QHry3rk9t1A86263cBtXPtHrwXrYxSDt71py+aDv0qtMUTMtQVpapO3HjuZxY+Uy/tWvnom1drVcNZfV94dx1JzTuO+UXBwtunrOgX4uDdnXfrm6H+jaoj3iVn2l9G1TXBxbu79320bZ26v3qVadpv65u626fadv6AF23QXV/kzQeBj9J0iL37XtkUVAqDzbPXdf5r6N+sFseZNfrtH1bXx4AVg8qTyquc4LHQllTUKrXgc5Bddm3Uj+jkvU6paagVD3YrS8fFo9INq3PonYaRlTKdsqAOchBdVmnn1G/MsRVP9N+wmy9nWqd6n5Qr1fvW3WUrGxrQTuVgFPvVxkO6tfDwWPbutxu9dMQ69uh3Edn1jcH7WqArNY9d90JzKxfWLYcxSwDcKk6al7dBvV+ddt3qiG0vANp2z5afuHStB+UdaojrOV+3fT3ULaz7qKFV+aVXxyU7dSDdnU/KNuqftlQX1bZVrewKGl4Bj9JUqMyKFVtXncqW9euaTzYLdXDVXmwWz2grD4Som0Er2yrutx6G03Bou3gvWyvrY1upy22jZLVl1nqtj7VUY7yYLfcXtUQUz2ohseCT/UguGkbdKvTOOpXfD71z7S6DeoH79V2ym29KPzV9p1qv/sZISpDT33fqqqGv3Lbldu0aVS2Gsrqy6rua3XlqZzVOvX9s3q659Ft0DByXv982oJSqb4+dU2jZPUb4TTt100he8vak1tHs9v2g/p2qLfbNJJb1qn/zVX/Hrqtj6ThGPwkSa0WjeBVRhJK5YF52+MGoH20q2nUr596sDgE1PtWjnQsGv1pOGWv16gf0PNgtzpiUV+f6ulq9XUq6zWdftftoLoaXuoHyFvXrukaFOq6bYPqOtX7DJ0Roqbrz8ptV63X1E59RLi6vctRsuo2qweW+qmoTcvstl7QfT/oFvLb6jR9EXLKppnGLzSa+lRfp+ry2v7Gmup02w/KOm397ta3UnW/7raftl0rWq9Tfh793shH0mAMfpKkrppGB9pOQYTHAklZp3qw23i7/Er5Bd/yV0aI2u5m2HZAWT8YbzpttHpw2Sso1U8rrYeR+vKqr7vVqZavr39T6CnVR/2q+gk9TafjVrdBfRSqWq+s2yvAlOvVpFvQhvbr1tr2g25fGjTto/WRtqa+Vet1+wy6bYNud6McNCj1+hKk2414yuV2q9O0vzVto6a/0259a7qusqmd+r7Sa30kDW65H+AuSTrGbV53KnP7H+150F8qD+h6PftvUEdPl6wcIG5edypbi9vs96pXrdPr8RRtB6Fb1p7ceU5gwzVITXW6jYr2aqdJ0zYoy9dvpDNonWrZ+qmpbeqnOjYpT/Frqrt53akLHtre9AD2+sPJoRNa6qev1vvSKxD3s+/UT51t+1uo37ETHtt2M+tPY0vx6IQ29XWsr09bMO/nMy3rlOtQffD9nj37R7qDZtlW/SyA+b0P9326Zrf1kTQaR/wkSV01fRtfBpGmQNJ2DVU/bTTdbbGt/NHXLe30G0yhdyhtGlFZcDrfgCOS1eXWT2/rZwSvafSm+rtNGR76+XyaHoXQ1E7PUb+GbVOt03RdHHTfdk3Xe1b3nV7PyKu3NYr6vjOOZ871+/n0umlNW71q+UH+5up1BhmVG6Rv1f5JGh9H/CRJraojJdXrofoZTel2UF0dcVgwvXK3wn6WVy5nkBt1tLVfPR20HurKNvodgaiO8pQHx9WRq6YbtcBj699Up77spv41qa/P0ZG1Lp9Pv59pdURpbv+jjSGo2lY/o6vw2HZo6hewaH37CQiD7qddT60cIDA2BcFBr0Gs9qk+qtbWRj/9KNvoNhrX+HdXq9PP+jSNtPfsW6XOurM29mxDUncGP0lST92CUpOmg+zytK9qGKhfx1c9oK8/f62+vLbl1E9vqweSfoJPXa9wVe9v0/tuy+23jbpe61Kue/3Aut/wVzfIdqu2XX529+17pPX0w7Z2y/ba1mUQbduhqteXBv2cNlqq7nttX3Z0C+bVfaPXurRN66atT219G6Qv1end9rduXxhIGi+DnySpq6YRsPq1av18Y9+0zHL+sEGk3xDXzyhTt35C9yDX72hS00hp4/V3LQf89RBVX2aT6oF1uR26BaimA/Gm/leX17SMbvV66da3fnUL/8P0rbqMQYPJoF829Bvo2953W27VoH0q60Dvv4dB9fpy4tsHXqKkusjMle5Dq9nZ2dy+fftKd0OSjjsf+MTfLprW76ll9YPQelCpL7PXsqq6HexW22kLCE3l2/rWptv69KMavnq1MUw7g4bcbne/7FZ+XHWa2q9vm7Z2duw/uGCErr6s+vy2PvXqWzmv7TNr6nvT5zyufaebbtuuqf9NAXaY/rXt171Cab9/D89++pMH7pOOLxFxW2bOrnQ/JpkjfpKkvvQzytF0vVfbKFHTqFPTgW314L0+v58RrHr5tn71q1qvn4PWphG3+utuhh2trN+dsrybZj0ItW3zurbPoFu46tZGt3lNbTTpPDy+eV5Zp61/w2zXXiGqabnj2OealjtMuaZgWx/JK99327fr8/od+e2330txwxxJjvhJkhqUI379BJuqpgO9YUaISt1urV8/mG9qp9e0fg0bvvpVrucgdyJtqr+SmvreFrr6CYtluapO0GvW9GzAnfuOHJ0+7LZtuyZxlJHSUtNoWL/Lawr2wxq2D72WN67lAPzMsy8eaVmafo749eaInySpp+UITd1HdRYf2NfDwyAjIt1Ga0YJqr30CmjjCnD1gNT2wPSlWkZ1PdrWqSn8NZVtCnvzexe+n1m/MORV65TTq8seJCSNYzSrn+X0uulNVdN26mffaVvvce/no4xSD/M3Lak/Bj9J0iLjOtht0s8BatPBfnVaeYA/yMF8/dEDbcrw1+/phv2El15GDWrdRsL6mT9qG/Vg1a9eI3/1NquBb+fu6LR9ZjK/97Hw17ac6jatt9vt1Nbl0u/+Nop+R9CrI9D1v7Fe/Rp0u5WfRfXU3GGWI6k3g58kaUn1OlDsFUrKg/2Z9YvrdDuYH8UgB+FtB8n96hVyh1UNSdVttxT168GqaT40B9peI7310T14LPSVr8vw16Y+Itit3er8bvtTry8Ceo1oVveXXiG0yThGddv61jatn751227dRn8HWY6k4Rj8JEl9GfSb+GFH9qqqB/NtAbDpYH5cN/OoLnNc5aD76YvDhLT6MpqmD7Lc+nL6qd8U/urr2TRqO0g/Wtsuwl/Pcj0Cat0w4aOf01wHKdekny8GmspU133QbTGIYUcrBw3nkgZj8JMk9dTvQWp1tKObbgeuvQ726yGk2+hf22mb1XJVvU7b7OdguV6mn4P0XgG3rU5Zpp+A1M9yuy2naV59WYOMVHYbBey2nOpoX316W/irbqu2dttG0CYlfPRzOm8/+2a3903XSFbndfvMJE0+g58kaWyGCXzlQXm/ozvVOk3Lrh+U9hv6ymlN4bV+s5C6+gHzuE/VLNWD3iDbrLrc6rYbZhn1uqOMUrYF5V6neDYuq8+Rv2o7/farn2X0E4j6GYmrThsmTA+ra+Du8TdQGnYbVOcZLKWlYfCTJC0yzlGOfkb3xhFgqu2VNxsZ9FqjpnnDnla3FEYJaUuxnOryRg1/41INh9UQOEwfu4WQfq7N7DWqOMiyjhW9RgX7/Xsy/EnjZ/CTJA2k7cBtkg7UmsJfv48KGGakZRTjDmErpS1Yta3fwCGsx2hfW52lCH/97htLtQ9VR8lHCdxLaZhTnrvVlzQ6g58kqS/DnBrXVmecI1e9Dnx73TWy3+nqrX7qZ69rBts+u0kLw/2MWi9XABtmlHzQO7z2cz3nIIb9m/JvURqvob9KiYitEfHxiLgzIm6PiNcW08+MiJsj4u7i9xnF9IiI/xIR90TE30XE08e1EpKkpTXMCMdyjpiVP219qRu0b03Lbys3LXbujqFG2eDY3F6j9mXQINbvNurVxjDLaLpra9u8bvV62bnviOFNmiCjjPgdAv5tZn42ItYCt0XEzcArgFsy8x0R8UbgjcAbgBcCFxY/lwDvKX5LkibMKAdrgzyiYamNcpoetB8g96rTNuJVHzWZpOBTNWzgG0Z9NGqpt0nTSGM1+Cz1Yy967VPLNQo6zL5dlhvH6Oa4RxUl9TZ08MvMeWC+eL0/Iu4ENgOXAc8uil0LfIJO8LsMeH9mJvDpiNgQETPFciRJx4GVOKgfNIgO8oiEfpc56PyVVA99g9wpc1Rt22XcQbSfx1YMGuD6qb9Uo6DLfdpsv2G13/rV6YY/aemM5Rq/iDgf+G7gVuDsMsxl5nxEPK4othnYUak2V0xbEPwi4mrgaoDzzjtvHN2TJE2A5Qo7gxw8juPUuUEs1bKXM5wtp3EEvmG3zaghZJjHlIxi0Oc6jrvt0rCBuT7fACiN38jBLyJOB/4Y+MXM3BfR+o9004xF/xJn5jZgG8Ds7Oz0/S8mSceZlRjdGuaUu0lQBp1+gkrTyFyprD/I8nq1NY3BspemEDLsTVWWwyTs0+MafXT0Txq/kYJfRKyhE/o+kJkfKiZ/vTyFMyJmgPuL6XPA1kr1LcDOUdqXJC2ffg7ajpVr2CaxX9Xg1i1o9TMKNuzpmt2W3fZ8vKWynNcZ9jKJ+8sk63Yd5aDLkTQ+Qwe/6AztXQPcmZm/VZl1A3Al8I7i90cq018TEdfRuanLXq/vk6TJNOwB1zgO1MqQUg0rw4w4tZ0yNvCdCZch8DSFnHEHn7bl1UcHB1nWINtjkM9wKdb9eByxXEkreeqppGajjPg9C3gZ8MWI+Hwx7c10At8HI+Iq4KvA5cW8G4EXAfcAjwCvHKFtSdKUqY94tU0bNgCO2qem94Oq930SRrVG6UNT3W7r2E8AW6ptYvhbfgY+abKMclfPT9J83R7ApQ3lE3j1sO1JkqbXMKNNTcZxYL+UYWwSgt5S67WOx8M2GEQ/4VmSxmEsd/WUJKmX5TjgX8nTCTX5JmHUb9BrNKv9Xer+Ny2/OtI+rpsFSVoZBj9J0iL9HJz2uv5upYJVr4NjA9/0eWD3Kjaeebivsr2udRxEfVlLsd/1e6OeQUNhW3/6mW4AlI5NBj9J0lCaDgQnRbcRk0n2wO5VAH2HmEHLT7NBwl+TfkNcP3c+Xep9rh7ymq6JXYl+SJpsBj9J0lQbxwHvqKGi3zb6ba9atp/y06y+LcZl3De9GbdhHuuxlP0YdsS0egppE4OlND4GP0mSuiiDxVKFq7bg0tRet5AzTAAaZH2qy1+uEFxtp5/tcTwH4JU27Kmm3gxIWj4GP0mSGjQFqfq0ejDppSzfb0hbqtGsUZe/lAGrPvI5aLuGv5XTbfTPACetPIOfJKmrcV5HNmqQ6Ra0RunfKAFoOdqZRE3bf5R9ZVLCsEbXdg2ipJVl8JMkLdLPaNdK6PdUx7bgMQnrMK26jdTVVUPiUvbFkb+VY9jTtIiIK4BfAt6fme8qpp0O/GfgecCjwIPA6zPz1hXraB8MfpKkqWPAm2x+PpKOIVcA3wN8ICJOz8xvAO8FvgRcmJlHIuKfAE9ayU72w+AnSZIkSc3K4esEIiK+DbgE+JnMPAKQmfcC965Q//p2wkp3QJIkaSk4sihpDD4EbAe2Z+Z+4DuBz2fmMXcuuSN+kiRJktQgM68Frl3pfoyDI36SJGkqeXMXSUvgduCpEXHM5ahjrsOSJOn41i3QbTzz8NEfSRq3zPxHOqd+/mpEBEBEXBgRl61sz3oz+EmSpGNGGeiawp1hT9IyeRVwDnBPRHwR+H1g58p2qTev8ZMkSUuuKZQ9sHtV12f6VeuUZZvKtM2TpKWQmfuAn1vpfgzK4CdJmmg+jH16VUfv+nno+zDzJEkdBj9J0iKjHEjXD+D7XdagozaD9tGg2K7cloNso7YRvH7LduuDQU6Sxs/gJ0kaq2EP2pf6YL/X8pc6GA46ctmtv73qNIXvpjqDXCNXrd9WbhyfoaFPkpaGwU+SJHqHnkECST8hqd/5w9RpuxauNMyomoFWZqqiAAAgAElEQVRMS+ncM/Po6527YwV7InWcserEPJRHhqr7jTx8U2a+YMxdGpnBT5KkHqYtJE16/3R8qYa+8r3hTyvtUB7haSevH6ruJ7+5e+OYuzMWBj9JkiQtq3rYa5q/c3f0LCepfwY/SZIkLYtBgpyhTxqvkYNfRKyi8/T6+zLzRyLiAuA64Ezgs8DLMvNARJwEvB94BvAg8JLM/PKo7UuSJB1PmgJR9dTIcv4gp0sOU6ffZUqaDOMY8XstcCewrnj/68A7M/O6iPg94CrgPcXvhzLzCRFxRVHuJWNoX5IkaWr1E6CaygwTvAxr0vQaKfhFxBbgh4G3A6+LiACeC7y0KHIt8Ct0gt9lxWuA64HfiYjITP+FkaQJM8rBX9PIQ3Xech9YepMITYpB71xpCJM0TqOO+L0L+CVgbfH+LGBPZh4q3s8Bm4vXm4EdAJl5KCL2FuUfqC4wIq4GrgY477zzRuyeJGm5dTtYXYkD2X7Cp+FwOqz0IwG8fk3SJBs6+EXEjwD3Z+ZtEfHscnJD0exj3mMTMrcB2wBmZ2f9V1GSNFbjOiVupZXBpq3v4xpRWulQXL21/7iCVbdR6fr8QZYraXqsjhM4e/VpQ9bePda+jMsoI37PAn40Il4EnEznGr93ARsiYnUx6rcF2FmUnwO2AnMRsRpYz6RuFUmSJlw/t8NfjnbGpdtpwOPuw3JtO0maJCcMWzEz35SZWzLzfOAK4GOZ+TPAx4EXF8WuBD5SvL6heE8x/2Ne3ydJksCwJUlLbejg18Ub6Nzo5R461/BdU0y/BjirmP464I1L0LYkSZIkTayIeF9E3B8Rf98y/2ci4u+Kn7+JiKeOo92xPMA9Mz8BfKJ4fS9wcUOZR4HLx9GeJGlpzawfrt783vG2PY7lSZI0Yf4r8Dt0nnHe5EvAD2TmQxHxQjr3P7lk1EbHEvwkSYLhA+NyLc8gKUlaaZn5VxFxfpf5f1N5+2k6900ZmcFPknTc6CdIluFwkLLSSmrbV7vtn93273Hu1/UR/G4j+uP+okeaElcBfzaOBRn8JEmqGOTgc9IPVOshtnqg3dR3g+zkGWUfG7buUu3X9eVO+t+Pjm9r4gTOOeHkYatvjIjtlffbikfWDSQinkMn+P2zYTtSZfCTJGlKDXqgvVwH4v0EzHpf6qNFg7a1nNeOGmik494DmTk7ygIi4inAe4EXZuaD4+iUwU+SJC2rYYLROEevDGaSJllEnAd8CHhZZv7vcS3X4CdJWuTcdcM/7WfnviMLllG+lyRJEBH/HXg2nVNC54C3AmsAMvP3gLfQeSze70YEwKFRRxDB4CdJGrN6aBwlRI6DwVOSNEky86d7zH8V8Kpxt2vwkyRNtUGCZzUk9qpnoJQkHUsMfpIkFQYJics9kmnQlCSNwuAnSdIxwFNmJUmjMPhJkqSexh0824LkuNqZ9qA67CnMg9Tdue8I5647ofEU6H6276CnS6/0lxtS1ZoMzsnpikrTtTaSpLHYunbN0HV37D84xp5oWi31Qb4h4jHDbouyXlP9cWxfPyNpeRn8JEljNUpoHBfDpyRJCxn8JElTZ5jwaViUJE0zg58kSYxvpHLH/oNdl2XAlCStBIOfJElj1CtALvWpsAZLSVITg58kSVNkqYKlgVKSjm0GP0mS1NO4A6VBUtIkW31C8riTh3wszL7x9mVcDH6SpEW2rD256/y5/Y8uU080rSbh7q9VBlFJ087gJ0kaWK9gOCqDpZbbpAXRSWMwlo59Bj9J0sRZ6mBZZ9CUujMYS8c+g58k6bg3TNA0LEqShhERLwB+G1gFvDcz31Gb/3jgfcAmYDfws5k5N2q7IwW/iNgAvBd4MpDAvwDuAv4IOB/4MvBTmflQRASdFXwR8Ajwisz87CjtS5K0UoYdlTQwStLxKyJWAe8Gng/MAZ+JiBsy845Ksd8E3p+Z10bEc4H/CLxs1LZHHfH7beDPM/PFEXEicCrwZuCWzHxHRLwReCPwBuCFwIXFzyXAe4rfkiQdN8Z5GqshUpKOORcD92TmvQARcR1wGVANfhcB/6Z4/XHgw+NoeOjgFxHrgO8HXgGQmQeAAxFxGfDsoti1wCfoBL/L6CTXBD4dERsiYiYz54fuvSRJxzFvsiNJE2ljRGyvvN+WmduK15uBHZV5cyweDPsC8JN0Btl+HFgbEWdl5oOjdGqUEb9/AuwC/t+IeCpwG/Ba4OwyzGXmfEQ8rijftJKbAYOfJEkTqFuwNBRKmmarV8OZZx0ervL9PJCZsy1zo2Fa1t7/O+B3IuIVwF8B9wGHhuvMY0YJfquBpwP/OjNvjYjfpnNaZ5t+VpKIuBq4GuC8884boXuSJGmpeMqqJA1lDthaeb8F2FktkJk7gZ8AiIjTgZ/MzL2jNjxK8JsD5jLz1uL99XSC39fLUzgjYga4v1K+60oCFMOg2wBmZ2cXBUNJkjRdlvvxHSUDp6QV8Bngwoi4gM5I3hXAS6sFImIjsDszjwBvonOHz5ENHfwy82sRsSMivj0z7wIupXNR4h3AlcA7it8fKarcALymuIDxEmCv1/dJ0mTavO7UsSznvn2PjGU50lIYV+AsA+SgyzN4SsefzDwUEa8BbqLzOIf3ZebtEfE2YHtm3kDnfin/MSKSzqmerx5H26Pe1fNfAx8o7uh5L/BK4ATggxFxFfBV4PKi7I10HuVwD53HObxyxLYlSRNumABZhsWyruFRk27YALkSI51z+x9d0K7hU1p+mXkjnWxUnfaWyuvr6ZxNOVYjBb/M/DzQdOHipQ1lkzGlVUnS9KqHxXGNPjYxVOp4Uw+bg4RPQ6J0bBt1xE+SpGNWU6isjzhWp0nHs5W6FlPSeBj8JEmqaAqDbaOOBkJJmk6rVgenn71muMp3jrcv42LwkyRpSP2ehmpAlCStNIOfJElLbJDrFA2JkqSlYPCTJGmCDHszGwOjJKkbg58kSVPAu59Kkrox+EmSpK7qodIgKEnHHoOfJEkayDhHFw2RkrQ8DH6SJGnFjCtEGiAljVOsWcXJZ5+x0t0YK4OfJEk65o0SIMvQ6CmtkqaZwU+SJB3X2kJj23QDoaRjkcFPkiRpAKM+cqNa3xApabkY/CRJkpZBU2AcZFSxW+A0QErHjoh4AfDbwCrgvZn5joYyPwX8CpDAFzLzpaO2a/CTJC0ys/60Be/n9z68Qj2Rjk9lyLtv3yN9jTA2jSI6sihNnohYBbwbeD4wB3wmIm7IzDsqZS4E3gQ8KzMfiojHjaNtg58kqaeZ9actCH/1YFhnUJTGY5jTSvsdWTQMSiviYuCezLwXICKuAy4D7qiU+Tng3Zn5EEBm3j+Ohg1+kqS+9Ap79bKGP2my9RolrDMoSn3bGBHbK++3Zea24vVmYEdl3hxwSa3+EwEi4v+jczror2Tmn4/aKYOfJGlJlEGx10ihAVFaecOeTipNqxPWrOHkc2aGrf5AZs62zIuGaVl7vxq4EHg2sAX464h4cmbuGbZD5UIlSVoyvUYKHR2Ujj2jPDdxEAZMTaE5YGvl/RZgZ0OZT2fmQeBLEXEXnSD4mVEaPmGUypIkjcMgp5FKOn5sXnfqsoVMaZl8BrgwIi6IiBOBK4AbamU+DDwHICI20jn1895RGzb4SZImQlP4m1l/2tEfSZKOdZl5CHgNcBNwJ/DBzLw9It4WET9aFLsJeDAi7gA+Drw+Mx8ctW1P9ZQkTYzqdYH1sNd0zaAkSceazLwRuLE27S2V1wm8rvgZG4OfJGnidBvhq85rC4H9lJEk6Xhi8JMkHbP6uXNovyHQsChJmmYjBb+I+DfAq+jcgvSLwCuBGeA64Ezgs8DLMvNARJwEvB94BvAg8JLM/PIo7UuSBP3fHGbYcgZBSTq+xOo1nLTpnJXuxlgNfXOXiNgM/AIwm5lPpvNwwSuAXwfemZkXAg8BVxVVrgIeyswnAO8sykmSNPG8uYwk6Vg36l09VwOnRMRq4FRgHngucH0x/1rgx4rXlxXvKeZfGhFNDzCUJGniGP4kSceyoYNfZt4H/CbwVTqBby9wG7CnuE0pdB4+uLl4vRnYUdQ9VJQ/q77ciLg6IrZHxPZdu3YN2z1JksbOx0tIko5VQ1/jFxFn0BnFuwDYA/wP4IUNRbOs0mXeYxMytwHbAGZnZxfNlyRpErSFP68HlCRNolFu7vI84EuZuQsgIj4E/FNgQ0SsLkb1tgA7i/JzwFZgrjg1dD2we4T2JUmaON4dVJI0iUa5xu+rwDMj4tTiWr1LgfLp8i8uylwJfKR4fUPxnmL+x4qHE0qSNJU8LVSSNCmGHvHLzFsj4no6j2w4BHyOzimafwpcFxG/Vky7pqhyDfAHEXEPnZG+K0bpuCRJx4rq8wYdEZQkrYSRnuOXmW8F3lqbfC9wcUPZR4HLR2lPkqRj2bAPl5ckLa9YvYYTN82sdDfGatTHOUiSpDHodcdQTxmVJI1ipBE/SZI0fv2EP0cIJUmDcMRPkqRjUK8RQG8sI0mTKSJeEBF3RcQ9EfHGhvk/HxFfjIjPR8QnI+KicbRr8JMk6RjVFu6q0wyAkjQ5ImIV8G46zz+/CPjphmD3h5n5XZn5NOA/Ab81jrYNfpIkHePqQa9XGUnSirkYuCcz783MA8B1wGXVApm5r/L2NGAsj8DzGj9JkqZAP8Gu+lgJSdKK2AzsqLyfAy6pF4qIVwOvA04EnjuOhg1+kiQdZwyAktTdCWvWcOLGc4etvjEitlfeb8vMbcXraCi/aEQvM98NvDsiXgr8n8CVw3amZPCTJOk4NbP+tEXhzzuHStLIHsjM2ZZ5c8DWyvstwM4uy7oOeM84OmXwkyRJPW8SYwiUpLH4DHBhRFwA3AdcAby0WiAiLszMu4u3PwzczRgY/CRJOo71e9MXTw+VpNFl5qGIeA1wE7AKeF9m3h4RbwO2Z+YNwGsi4nnAQeAhxnCaJxj8JEnSAJpOD5Uk9S8zbwRurE17S+X1a5eiXYOfJEkaSD+jhIZDSZosBj9JkjR2beHQQChJK8MHuEuSpGXjg+QlaWU44idJkpaV1wlKmnSx+kRO2jT0c/wm0kQHvy9/+cvMzrY9AkOStFQe/da3VroLknTUySedtNJd0OR7+kp3YNJNdPA7//zz2b59e++CkqSxuusfv7TSXZCko7792y5Y6S5owkXEZ1e6D5POa/wkSZIkacoZ/CRJkiRpyhn8JEmSJGnKGfwkSZIkacpN9M1dJEkrY2bD2tZ583v2L2NPJElafrF6jY9zkCQd37qFwuNJUwCe2bDWYCxJmkgGP0mShtAWgA3GkqRJ5DV+kiRJkjTlDH6SJEmStEwi4gURcVdE3BMRb2yYf1JE/FEx/9aIOH8c7Rr8JEmSJGkZRMQq4N3AC4GLgJ+OiItqxa4CHsrMJwDvBH59HG0b/CRJkiRpeVwM3JOZ92bmAeA64LJamcuAa4vX1wOXRkSM2rDBT5IkSZLGZ2NEbK/8XF2ZtxnYUXk/V0yjqUxmHgL2AmeN2inv6ilJkiRJFSesXsMpm2aGrf5AZs62zGsaucshygzMET9JkiRJWh5zwNbK+y3AzrYyEbEaWA/sHrVhg58kSZIkLY/PABdGxAURcSJwBXBDrcwNwJXF6xcDH8vMkUf8PNVTkiRJkpZBZh6KiNcANwGrgPdl5u0R8TZge2beAFwD/EFE3ENnpO+KcbRt8JMkSZKkZZKZNwI31qa9pfL6UeDycbfrqZ6SJEmSNOUMfpIkSZI05TzVU5IkSZIqDh4+zPye/SvdjbEy+EmSJK2wb+6aX+kuTLR1Z21c6S5IxzyDnyRpRfU64D1l04wHxZIkjcjgJ0laZJKC1iT1RZKkY5U3d5EkSZKkKWfwkyRJkqQp1zP4RcTWiPh4RNwZEbdHxGuL6WdGxM0RcXfx+4xiekTEf4mIeyLi7yLi6ZVlXVmUvzsirly61ZIkSZIklfoZ8TsE/NvMfBLwTODVEXER8Ebglsy8ELileA/wQuDC4udq4D3QCYrAW4FLgIuBt5ZhUZIkSZK0dHre3CUz54H54vX+iLgT2AxcBjy7KHYt8AngDcX092dmAp+OiA0RMVOUvTkzdwNExM3AC4D/Psb1kSRJkqSRHDx8hPm9D690N8ZqoGv8IuJ84LuBW4Gzi1BYhsPHFcU2Azsq1eaKaW3T621cHRHbI2L7rl27BumeJEmSJKlB38EvIk4H/hj4xczc161ow7TsMn3hhMxtmTmbmbObNm3qt3uSJEmSpBZ9Bb+IWEMn9H0gMz9UTP56cQonxe/7i+lzwNZK9S3Azi7TJUmSJOm413YDzVqZx0fEbRHx+eLmmz/fz7L7uatnANcAd2bmb1Vm3QCUd+a8EvhIZfrLi7t7PhPYW5wKehPwgxFxRrECP1hMkyRJkiS130Czah74p5n5NDo3znxjRJzba8E9b+4CPAt4GfDFiPh8Me3NwDuAD0bEVcBXgcuLeTcCLwLuAR4BXgmQmbsj4t8DnynKva280YskSZIkqfUGmkdl5oHK25Po8yzOfu7q+Umar88DuLShfAKvblnW+4D39dMxSZIkSToGbYyI7ZX32zJzW591F9xAMyIe11QoIrYCfwo8AXh9Zva8hK6fET9JkiRJOm4cPHyE+/Y9Mmz1BzJztm1mRHwUOKdh1i/320Bm7gCeUpzi+eGIuD4zv96tjsFPkiRJkpZJZj6vbV5EfD0iZorRvuoNNNuWtTMibge+D7i+W9mBnuMnSZIkSVoybTfQPCoitkTEKcXrM+jck+WuXgs2+EmSJEnSZHgH8PyIuBt4fvGeiJiNiPcWZZ4E3BoRXwD+EvjNzPxirwV7qqckSZIkTYDMfJDmG2huB15VvL4ZeMqgy3bET5IkSZKmnMFPkiRJkqacp3pKkiRJUsWBI8nc/kdXuhtj5YifJEmSJE05g58kSZIkTTmDnyRJkiRNOYOfJEmSJE05g58kSZIkTTmDnyRJkiRNOYOfJEmSJE05n+MnSZIkSRUHDic79h9c6W6MlSN+kiRJkjTlDH6SJEmSNOUMfpIkSZI0ASLizIi4OSLuLn6f0VLuvIj4i4i4MyLuiIjzey3b4CdJkiRJk+GNwC2ZeSFwS/G+yfuB38jMJwEXA/f3WrDBT5IkSZImw2XAtcXra4EfqxeIiIuA1Zl5M0BmfiMzH+m1YIOfJEmSJE2GszNzHqD4/biGMk8E9kTEhyLicxHxGxGxqteCez7OISJOBv4KOKkof31mvjUiLgCuA84EPgu8LDMPRMRJdIYenwE8CLwkM79cLOtNwFXAYeAXMvOmXu1LkiRJ0nI6eDjZue/IsNU3RsT2yvttmbmtfBMRHwXOaaj3y30ufzXwfcB3A18F/gh4BXBNr0q9fAt4bmZ+IyLWAJ+MiD8DXge8MzOvi4jfoxPo3lP8figznxARVwC/DrykGJK8AvhO4FzgoxHxxMw83OcKSpIkSdKkeyAzZ9tmZubz2uZFxNcjYiYz5yNihuZr9+aAz2XmvUWdDwPPpEfw63mqZ3Z8o3i7pvhJ4LnA9cX06vmn1fNSrwcujYgopl+Xmd/KzC8B99C5EFGSJEmSBDcAVxavrwQ+0lDmM8AZEbGpeP9c4I5eC+7rGr+IWBURn6eTOG8G/hHYk5mHiiJzwObi9WZgB0Axfy9wVnV6Qx1JkiRJOt69A3h+RNwNPL94T0TMRsR7AYozJv8dcEtEfBEI4Pd7LbifUz3LhT8tIjYAfwI8qalY8Tta5rVNXyAirgauBjjvvPP66Z4kSZIkHfMy80Hg0obp24FXVd7fDDxlkGUPdFfPzNwDfILOOaQbIqIMjluAncXrOWArQDF/PbC7Or2hTrWNbZk5m5mzmzZtqs+WJEmSJA2oZ/CLiE3FSB8RcQrwPOBO4OPAi4ti1fNPq+elvhj4WGZmMf2KiDipuCPohcDfjmtFJEmSJEnN+jnVcwa4tng2xAnABzPzf0XEHcB1EfFrwOd47C4y1wB/EBH30BnpuwIgM2+PiA/SufDwEPBq7+gpSZIkSUuvZ/DLzL+j84yI+vR7abgrZ2Y+Clzesqy3A28fvJuSJEmStDwOHob5vSvdi/Ea6Bo/SZIkSdKxx+AnSZIkSVOur8c5SJKOL6dsmhnLcr65a36gdnqVlyRJwzH4SZKWzKABclyBc7lUg2pT3/sNsuNY72M9NFe3QT/r4pcGkjQYg58kSUPqFdiWM8gea6G5m2HWZZrWX5KWgtf4SZIkSdKUc8RPkiRJkioOHoKdu2OluzFWjvhJkiRJ0pQz+EmSJEnSlDP4SZIkSdKUM/hJkiRJ0gSIiDMj4uaIuLv4fUZLuV+PiL8vfl7Sz7INfpIkSZI0Gd4I3JKZFwK3FO8XiIgfBp4OPA24BHh9RKzrtWCDnyRJkiRNhsuAa4vX1wI/1lDmIuAvM/NQZj4MfAF4Qa8FG/wkSZIkaXw2RsT2ys/VA9Q9OzPnAYrfj2so8wXghRFxakRsBJ4DbO21YJ/jJ0mSJEkVhw4HD+xeNWz1BzJztm1mRHwUOKdh1i/3s/DM/IuI+B7gb4BdwKeAQ73qGfwkSZIkaZlk5vPa5kXE1yNiJjPnI2IGuL9lGW8H3l7U+UPg7l7teqqnJEmSJE2GG4Ari9dXAh+pF4iIVRFxVvH6KcBTgL/otWBH/CRJkiRpMrwD+GBEXAV8FbgcICJmgZ/PzFcBa4C/jgiAfcDPZqanekqSJEnSsSAzHwQubZi+HXhV8fpROnf2HIinekqSJEnSlDP4SZIkSdKU81RPSZIkSapYdSg584Gel801emTMfRmXvkf8irvHfC4i/lfx/oKIuDUi7o6IP4qIE4vpJxXv7ynmn19ZxpuK6XdFxA+Ne2UkSZIkSYsNcqrna4E7K+9/HXhnZl4IPARcVUy/CngoM58AvLMoR0RcBFwBfCfwAuB3I2LopyJKkiRJkvrTV/CLiC3ADwPvLd4H8Fzg+qLItcCPFa8vK95TzL+0KH8ZcF1mfiszvwTcA1w8jpWQJEmSJLXrd8TvXcAvAUeK92cBeyrPi5gDNhevNwM7AIr5e4vyR6c31JEkSZIkLZGewS8ifgS4PzNvq05uKJo95nWrU23v6ojYHhHbd+3a1at7kiRJkqQe+hnxexbwoxHxZeA6Oqd4vgvYEBHlXUG3ADuL13PAVoBi/npgd3V6Q52jMnNbZs5m5uymTZsGXiFJkiRJ0kI9g19mvikzt2Tm+XRuzvKxzPwZ4OPAi4tiVwIfKV7fULynmP+xzMxi+hXFXT8vAC4E/nZsayJJkiRJajTKc/zeAFwXEb8GfA64pph+DfAHEXEPnZG+KwAy8/aI+CBwB3AIeHVmHh6hfUmSJEkau1WHDrB2z9xKd2OsBgp+mfkJ4BPF63tpuCtnZj4KXN5S/+3A2wftpCRJkiRpeIM8x0+SJEmSdAwy+EmSJEnSBIiIyyPi9og4EhGzLWW2RsTHI+LOouxr+1m2wU+SJEmSJsPfAz8B/FWXMoeAf5uZTwKeCbw6Ii7qteBRbu4iSZIkSRqTzLwTIKLpEehHy8wD88Xr/RFxJ7CZzk00WzniJ0mSJEnHoIg4H/hu4NZeZR3xkyRJkqSKPHKAI9/40rDVN0bE9sr7bZm5rXwTER8Fzmmo98uZ+ZGG6Y0i4nTgj4FfzMx9vcob/CRJkiRpfB7IzMYbswBk5vNGbSAi1tAJfR/IzA/1U8dTPSVJkiTpGBGdCwCvAe7MzN/qt57BT5IkSZImQET8eETMAd8L/GlE3FRMPzcibiyKPQt4GfDciPh88fOiXsv2VE9JkiRJmgCZ+SfAnzRM3wm8qHj9SaD9tp8tHPGTJEmSpCln8JMkSZKkKeepnpIkSZJUkUcOcODhr650N8bKET9JkiRJmnIGP0mSJEmacgY/SZIkSZpyBj9JkiRJmnIGP0mSJEmacgY/SZIkSZpyBj9JkiRJmnI+x0+SJEmSKo4cPsA398+tdDfGyhE/SZIkSZpyBj9JkiRJmnJ9Bb+I+HJEfDEiPh8R24tpZ0bEzRFxd/H7jGJ6RMR/iYh7IuLvIuLpleVcWZS/OyKuXJpVkiRJkqRjT0RcHhG3R8SRiJjtUm5RPutlkBG/52Tm0zKz7MAbgVsy80LgluI9wAuBC4ufq4H3FJ07E3grcAlwMfDWMixKkiRJkvh74CeAv+qjbD2fdTXKqZ6XAdcWr68Ffqwy/f3Z8WlgQ0TMAD8E3JyZuzPzIeBm4AUjtC9JkiRJUyMz78zMu5Zi2f0GvwT+IiJui4iri2lnZ+Z80cF54HHF9M3AjkrduWJa2/QFIuLqiNgeEdt37drV/5pIkiRJ0vGhKZ911e/jHJ6VmTsj4nHAzRHxD13KRkvH2qYvnJC5DdgGMDs7u2i+JEmSJC2lQwf33/S1uVs2Dln95Np1d9uKjANARHwUOKeh3i9n5kf6bGNRPsvMrqeH9hX8MnNn8fv+iPgTOtfofT0iZjJzvjiV8/6i+BywtVJ9C7CzmP7s2vRP9NO+JEmSJC2XzFyyS9Iy83ljWEZTPusa/Hqe6hkRp0XE2vI18IN0Ljq8ASjvzHklUKbTG4CXF3f3fCawtzgV9CbgByPijOKmLj9YTJMkSZIk9aFLPuuqnxG/s4E/iYiy/B9m5p9HxGeAD0bEVcBXgcuL8jcCLwLuAR4BXgmQmbsj4t8DnynKvS0zd/e5fpIkSZI01SLix4H/G9gE/GlEfD4zfygizgXem5kvoiWf9Vx25uReRjc7O5vbt/f1WApJ0hjte/CBle6CJB217qxhL7XS8SIibuv3sQbHq1Ee5yBJkiRJOgYY/CRJkiRpyhn8JEmSJGnKGfwkSZIkacoZ/CRJkiRpyhn8JEmSJGnKGfwkSZIkacoZ/CRJkiRpyhn8JEmSJGnKGfwkSZIkacoZ/CRJkiRpyhn8JEmSJGnKGfwkSZIkacoZ/CRJkiRpyhn8JEmSJGnKGfwkSZIkacpFZrrYwckAAAvMSURBVK50H1pFxH7grpXuhybCRuCBle6EVpz7gcD9QB3uBwL3Az3m8Zm5aaU7MclWr3QHergrM2dXuhNaeRGx3X1B7gcC9wN1uB8I3A+kQXiqpyRJkiRNOYOfJEmSJE25SQ9+21a6A5oY7gsC9wN1uB8I3A/U4X4g9Wmib+4iSZIkSRrdpI/4SZIkSZJGZPCTJEmSpCln8JMkSdLEiog1K90HaRoY/CRJEykiYqX7IGnlRMQzI+I64Dci4skr3R/pWLdiwS8inhYRPxcR56xUHzQZ3BcE7gfqiIjviogXR8Qp6d3Hjlv+e6CIuBx4D/C/gJOB1xXT/UJIGtKy39WzGK7/HWAWuBP4FrAtM29d1o5oxbkvCNwP1BERJ9HZD74H+ApwL/DOzPzqinZMy8p/D1SKiDcDJ2fmWyLi24E3AP9HZh5c4a5Jx6yVGPH7LmB9Zj4jM3+26MMDK9APrTz3BQF8J+4Hgh+gsx88DfgXwBOBR1a2S1oBT8Z/D45LEfFTEfG6iPjeYtJdwE9ExC8BnwLOBd4dEd+zYp2UjnHLEvwi4vkR8fzi7WHgpyJifUT8BPBM4NKI+O6irEP4U6w4hetfFW/dF45TEfETEfGu4u0JuB8clyLi6cU3+QAHgOcUr58NrAeeGxFbVqJvWj4RcUFEnFy+xX8PjisRsSoi3kJnRA/g9yPiR4EPAa8Fvh94eWa+ANgF/KSnAEvDWdJTPf//9u499uq6juP48yVoaGRempcEi+YlrVRMCXVmKqa4VfOOaxPNTDJTi0xFyYQmJmjKXMu7ztSZl9Kcc4VZis6ad4RMhy4vqcyYhpgI8uqPz+cUMpr8fr/z4/x+39/rsTE453zP4XP2fX8/3/P+XCV9BpgMbAqcavvh+vy5wCeA/YCJlJbdnYCJtp/ptQJFx0gaClwLbAFcDNxse5mkacCWJBYGBEnbA2cB2wE7AsNs/0PSecBwEgcDgqQRlOF8G1ES/0m275E0k3Jv2A2YREkA3wDOs/1Sh4obvUTSJylzuNahnOfJtudJmgpsBexL6oMBQdL1wOW2/yjpEOBEYILtv0m6BTjD9rOSRgE/Bo6wvaiDRY7ol9re49dqjZO0EXAfsND23rYfltT6/84EngYOtX0dcBHwPLBHu8sTnbNSy+xw4DXbo23fCLRaHCZR5nEkFhpqhTrhi8DlwEO2R1LOdWtIzxmUODgkcdBMK9UHPwAet70bcDvwzfr89yjn/cu2rwCmAR8CtiUaYRVx8Gfb+wL3AlMlbWN7MjCX3BcaS9JRkvaStEF96jVgQ0mDbd9KOf/jas/efODQetxI4J01X+KIZuiNoZ5DAGwvBKZTbtpIOhoYI2nrulLbMuDweuw/KT1B83qhPNE5Q1b49w7AMIA61PMsSfvUWHgLGAeJhYZat/49j/KDfqakdYCtgaUANQ7eI3VCkw2B//7wX0w998D6wDxJ29t+jzKf6wAA23MpjUbp7WuOVhwMro/nAti+BPg88I3aSPwOqQ8aRcXmku4FxgNfp8zZG0q57j8HDK2HXwIcRLkv/A7YVdJDwGGUEQLp7YvohrYlfnUe3+8pe62Mq09fTLlYXwG+ChwI3CZpK+AOYKyk6ZLup/wIeK5d5YnOWSEWzpd0ZH36UeAVSVdRenneBE6TNIHS0jtG0ozEQnOsFAfjbL9ue7GkIbbfBeZQbvwtdwAHpk5olpXuDYfXJH82sLWkxyhJ3iDgGkljKYnAIZKm1DhYACzI3K7+bRVxsAxYCIyUtKOkHYGnKIn+NsBdpD5oDEmD6rX/EeDl2st7AuW3wMXAzyk9ujtIWs/208CzwGG27wGOAo6zPaa+FhHd0JY5fjWRux44F3iBMnxjru1zJX0F2Mj2tfXYqykX/Vl1svYoYIHtX/e4INFx/ycWHqVU7D+lzNkZbXuppKMoE/dPoqzsOJrEQiOsIg4mAvNqnbB2Pf97AccB37e9oL5vB0rDQOKgAVYRB6cCj9qeURd1mWb74Hrs2cBQ26dK2hPYB5hj+7YOFT/aZBVx8EPgfuA64BTgC8AG9d8TgQdsXyJpp/pa6oN+qvbsTqE07txF6eE/1Pb4+vpawKuUhZ12oZzvP9m+qc77m5mtPCLaZ/AHH7Jqrfl6tpdTLtRHbN9eX7sHuFDSFbZ/2zq+HnsnMKa+9zHgsZ59hei0D4iFPwAXAFdS5vKMpAzVuAF4AjgYWG77ifo4+qnViINWnbCgvmVtysqNC1ufYftJ4Mk1WvBoqw+Ig1mUOLiOct5flLSd7b8C9wCn1HvF/ZTEIPqp1YiDCyiLfE2V9Cnbz9XXHqDO4bL9OPB4J8ofPVcb92YCD1IagKdSkv+9JY2y/RfbyyWdA0y3faCkN4HjJJ0OPEMZGRIRbdKtoZ6SjqHMuZhan5oDHFlX6ILyg24+MKP1nnpxjwfOBu7uZnmjj1mNWBhMmZR/vu37KJP0J0o6DbiRMuQrS3T3c92sE2ZRWnh3X2MFjV61mnHwXH19EWVVz5MknQxcCswCnPqgf1vN+8J84Gf18fP1fd+i7OH46Joqa/Sq5cAM29+2fTllKO8I4EeU1VxbDQS3Am9LGm77N8CxlIW+jrCdvTwj2qjLiV+dhPs1yrC9sZI+XVvprwWm1da6PYGjgY0lbSJpY0nTgWOAY1utftG/dSEWjgKGSdrM9h3ABMrS3cfanmF7uXtzX5HoVd2oEzar71ub0hD0QkcKHm3VhTgYT5nHNYjS+j8P2JlSH/zCVUe+RPRYF+uDDSVtatuSTqEM/T7edhK/ZngE+JWkQfXxA8CWtq8BBkn6bu0RHgYstf0igO1XWz3AEdFe3ZrjJ2lL2y+o7L01wvYR9cL+KLC97dmShlNa+1rLdG9h++9tK3n0CV2MhQm2swxzA3UhDqZQ4mBJRwscvaILcfATykIN73a0wNErunhfON72krqgR3p3GkzSNcCTti+sazwcR9nHd1vgMtvTO1m+iIGgW0M9bbda6C8CRkjavy7D/abt2fW1CcDb9fhlSfqaqYuxsHRVnxH9Xxfi4N+UrVyigboQB4spy7RHA3XxvrCsvidJX0NJGlSHdG5KWb0ZylDvScB5wJeS9EWsGT3azsH2q5RFOybVx+9JGiWptYjHFJclm6PhVjMW8kOv4RIHAYmDKBIHUS2nzO99nbJdw53AZMrCbrNtv9zR0kUMID3azqG1UqekW4BXgCWUyfnP2p7fpjJGP5BYCEgcRJE4CEgcxP9IGk1Z3fNB4GrbV3a4SBEDUk97/JZLWg/YBDgSeMH23anQB57EQkDiIIrEQUDiIN7nJeBMYN8kfRGd0+19/FZwAmXp5f2yYMOAl1gISBxEkTgISBwEYPslYFqnyxEx0PVoqCe8b2P2GOASCwGJgygSBwGJg4iIvqTHiV9ERERERET0bT2a4xcRERERERF9XxK/iIiIiIiIhkviFxERERER0XBJ/CIiouNUzJY0doXnDpd0dyfLFRER0RRZ3CUiIvoESZ8FbgZGAoOAx4EDerLvm6TBtpe1qYgRERH9VhK/iIjoMySdDywGPgwssj1V0njgO8A6wIPAiXVz8MuAnYF1gZtsT6mf8RJwKXAAcJHtmzvwVSIiIvqUdmzgHhER0S7nUDb8fhfYpfYCHgTsbntZTfbGATcAp9teKGkwcK+kW2zPq5+z2PYenfgCERERfVESv4iI6DNsL5Z0E/CW7SWSxgC7Ag9LgtK792I9/EhJx1LuZR8Htgdaid9Na7bkERERfVsSv4iI6GuW1z8AAq6yPXnFAyRtDZwMjLL9hqRfAkNWOGTxGilpREREP5FVPSMioi+bBRwu6WMAkjaWtCWwPrAI+JekzYH9O1jGiIiIPi89fhER0WfZniPpHGCWpLWApcAE4GHKsM6ngOeABzpXyoiIiL4vq3pGREREREQ0XIZ6RkRERERENFwSv4iIiIiIiIZL4hcREREREdFwSfwiIiIiIiIaLolfREREREREwyXxi4iIiIiIaLgkfhEREREREQ2XxC8iIiIiIqLh/gN+Io+CHgWpdgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x432 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots(2,1,sharex=True,figsize=(14,6))\n",
    "plt.subplots_adjust(hspace=0)\n",
    "levs = np.arange(-1.5,1.5,0.1)\n",
    "\n",
    "# 01deg\n",
    "t_up = temp_hov_01deg.sel(st_ocean=slice(0,1000))\n",
    "t_lo = temp_hov_01deg.sel(st_ocean=slice(1000,6000))\n",
    "\n",
    "axx=ax[0]\n",
    "t_up.T.plot.contourf(ax=axx,levels=levs,cmap = cm.cm.balance,yincrease=False,add_colorbar=False)\n",
    "axx.set_ylabel('')\n",
    "axx.set_xlabel('')\n",
    "plt.xlim([pd.datetime(1958,1,1),pd.datetime(2017,12,31)])\n",
    "axx.set_title('(c) ACCESS-OM2-01')\n",
    "axx=ax[1]\n",
    "p1=t_lo.T.plot.contourf(ax=axx,levels=levs,cmap = cm.cm.balance,yincrease=False,add_colorbar=False)\n",
    "axx.set_ylabel('')\n",
    "axx.set_xlabel('Year')\n",
    "plt.xlim([pd.datetime(1958,1,1),pd.datetime(2017,12,31)])\n",
    "\n",
    "ax1 = plt.axes([0.92,0.25,0.015,0.5])\n",
    "cb = plt.colorbar(p1,cax=ax1,orientation='vertical')\n",
    "ax1.xaxis.set_label_position(\"top\")\n",
    "cb.ax.set_xlabel('°C')\n",
    "\n",
    "#savefigure('tz_hovmoller')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:analysis3-18.10]",
   "language": "python",
   "name": "conda-env-analysis3-18.10-py"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}