override_variable_names_in_xorca_v0.2.0.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preamble"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/sfs/fs6/home-geomar/smomw122/miniconda3_20180820/envs/xorca_env/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n"
     ]
    }
   ],
   "source": [
    "from xorca.lib import load_xorca_dataset\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get all file names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = Path(\"original_data/\")\n",
    "data_files = list(sorted(data_path.glob(\"a13c*.nc\")))\n",
    "aux_files = list(sorted(data_path.glob(\"[m, n]*nc\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([PosixPath('original_data/a13c_1m_20000101_20001231_grid_U.nc'),\n",
       "  PosixPath('original_data/a13c_1m_20000101_20001231_grid_V.nc'),\n",
       "  PosixPath('original_data/a13c_1m_20010101_20011231_grid_U.nc'),\n",
       "  PosixPath('original_data/a13c_1m_20010101_20011231_grid_V.nc')],\n",
       " [PosixPath('original_data/mesh_mask_nemo.Ec3.2_O1L75.nc')])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_files, aux_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load dataset\n",
    "\n",
    "We currently use names according to the standard where what you call `\"uo\"` is `\"vozocrtx\"` etc.  All variable names are stored in dicts  in `xorca.orca_names`.  To (partially) override a dictionary from `xorca.orca_names`, add a keyword arg to `load_xorca_dataset()` called `update_<dict_name>`.  We add variables `uo`, `u2o`, `umo`, `tauuo` living on the U grid `(y_c, x_r)` and the analog for the V grid `(y_r, x_c)`.\n",
    "\n",
    "We also choose chunks that are fitting to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_xorca = load_xorca_dataset(data_files=data_files, aux_files=aux_files,\n",
    "                              decode_cf=True, \n",
    "                              input_ds_chunks={\"time_counter\": 4, \"t\": 4,\n",
    "                                               \"deptht\": 25, \"depthu\": 25,\n",
    "                                               \"depthv\": 25, \"depthw\": 25,\n",
    "                                               \"y\": 100, \"x\": 60},\n",
    "                              target_ds_chunks={\"t\": 4,\n",
    "                                                \"z_c\": 25, \"z_l\": 25,\n",
    "                                                \"y_c\": 290, \"y_r\": 290,\n",
    "                                                \"x_c\": 60, \"x_r\": 60},\n",
    "                              update_orca_variables={\n",
    "                                  'uo': {'dims': ['t', 'z_c', 'y_c', 'x_r']},\n",
    "                                  'u2o': {'dims': ['t', 'z_c', 'y_c', 'x_r']},\n",
    "                                  'umo': {'dims': ['t', 'z_c', 'y_c', 'x_r']},\n",
    "                                  'vo': {'dims': ['t', 'z_c', 'y_r', 'x_c']},\n",
    "                                  'v2o': {'dims': ['t', 'z_c', 'y_r', 'x_c']},\n",
    "                                  'vmo': {'dims': ['t', 'z_c', 'y_r', 'x_c']},\n",
    "                                  'tauuo': {'dims': ['t', 'y_c', 'x_r']},\n",
    "                                  'tauvo': {'dims': ['t', 'y_r', 'x_c']}})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (t: 24, x_c: 360, x_r: 360, y_c: 290, y_r: 290, z_c: 75, z_l: 75)\n",
       "Coordinates:\n",
       "  * z_c      (z_c) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ...\n",
       "  * z_l      (z_l) float64 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 ...\n",
       "  * y_c      (y_c) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ...\n",
       "  * y_r      (y_r) float64 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 ...\n",
       "  * x_c      (x_c) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ...\n",
       "  * x_r      (x_r) float64 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 ...\n",
       "    depth_c  (z_c) float64 dask.array<shape=(75,), chunksize=(25,)>\n",
       "    depth_l  (z_l) float64 dask.array<shape=(75,), chunksize=(25,)>\n",
       "    llat_cc  (y_c, x_c) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    llat_cr  (y_c, x_r) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    llat_rc  (y_r, x_c) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    llat_rr  (y_r, x_r) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    llon_cc  (y_c, x_c) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    llon_cr  (y_c, x_r) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    llon_rc  (y_r, x_c) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    llon_rr  (y_r, x_r) float32 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "  * t        (t) datetime64[ns] 2000-01-16T12:00:00 2000-02-15T12:00:00 ...\n",
       "    e1t      (y_c, x_c) float64 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    e2t      (y_c, x_c) float64 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    e3t      (z_c, y_c, x_c) float64 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "    e1u      (y_c, x_r) float64 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    e2u      (y_c, x_r) float64 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    e3u      (z_c, y_c, x_r) float64 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "    e1v      (y_r, x_c) float64 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    e2v      (y_r, x_c) float64 dask.array<shape=(290, 360), chunksize=(290, 60)>\n",
       "    e3v      (z_c, y_r, x_c) float64 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "    e3w      (z_l, y_c, x_c) float64 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "    tmask    (z_c, y_c, x_c) int8 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "    umask    (z_c, y_c, x_r) int8 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "    vmask    (z_c, y_r, x_c) int8 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "    fmask    (z_c, y_r, x_r) int8 dask.array<shape=(75, 290, 360), chunksize=(25, 290, 60)>\n",
       "Data variables:\n",
       "    uo       (t, z_c, y_c, x_r) float32 dask.array<shape=(24, 75, 290, 360), chunksize=(4, 25, 290, 60)>\n",
       "    u2o      (t, z_c, y_c, x_r) float32 dask.array<shape=(24, 75, 290, 360), chunksize=(4, 25, 290, 60)>\n",
       "    umo      (t, z_c, y_c, x_r) float32 dask.array<shape=(24, 75, 290, 360), chunksize=(4, 25, 290, 60)>\n",
       "    tauuo    (t, y_c, x_r) float32 dask.array<shape=(24, 290, 360), chunksize=(4, 290, 60)>\n",
       "    vo       (t, z_c, y_r, x_c) float32 dask.array<shape=(24, 75, 290, 360), chunksize=(4, 25, 290, 60)>\n",
       "    v2o      (t, z_c, y_r, x_c) float32 dask.array<shape=(24, 75, 290, 360), chunksize=(4, 25, 290, 60)>\n",
       "    vmo      (t, z_c, y_r, x_c) float32 dask.array<shape=(24, 75, 290, 360), chunksize=(4, 25, 290, 60)>\n",
       "    tauvo    (t, y_r, x_c) float32 dask.array<shape=(24, 290, 360), chunksize=(4, 290, 60)>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds_xorca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate the barotropic stream function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rename vars?\n",
    "\n",
    "As the var names differ from our standard, let's check what to rename by looking at the source code of `calculate_psi`. Load and then use the Ipython `??` notation to get a full listing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xorca.calc import calculate_psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mcalculate_psi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mSource:\u001b[0m   \n",
       "\u001b[0;32mdef\u001b[0m \u001b[0mcalculate_psi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m\"\"\"Calculate the barotropic stream function.\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m    Parameters\u001b[0m\n",
       "\u001b[0;34m    ----------\u001b[0m\n",
       "\u001b[0;34m    ds : xarray dataset\u001b[0m\n",
       "\u001b[0;34m        A grid-aware dataset as produced by `xorca.lib.preprocess_orca`.\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m    Returns\u001b[0m\n",
       "\u001b[0;34m    -------\u001b[0m\n",
       "\u001b[0;34m    psi : xarray data array\u001b[0m\n",
       "\u001b[0;34m        A grid-aware data array with the barotropic stream function in `[Sv]`.\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m    \"\"\"\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mgrid\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mxgcm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGrid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mperiodic\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"Y\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"X\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mU_bt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvozocrtx\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0me3u\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"z_c\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mpsi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgrid\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcumsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m \u001b[0mU_bt\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0me2u\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"Y\"\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m1.0e6\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mpsi\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0mpsi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my_r\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx_r\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# normalize upper right corner\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mpsi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpsi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrename\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"psi\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0mpsi\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mFile:\u001b[0m      ~/miniconda3_20180820/envs/xorca_env/lib/python3.6/site-packages/xorca/calc.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "calculate_psi??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rename `uo` to `vozocrtx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_xorca = ds_xorca.rename({\"uo\": \"vozocrtx\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 932 ms, sys: 286 ms, total: 1.22 s\n",
      "Wall time: 1.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "psi = calculate_psi(ds_xorca)\n",
    "psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 42.1 s, sys: 1min 37s, total: 2min 19s\n",
      "Wall time: 10.7 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAG5CAYAAAC0v4EiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xv8HVdd7//3JxTbBrAtTYFeCZKCgRAqxJIjR0FatHKUouI5oD8oFa16QKp44/JTih74IT+EH96QKG1BgVoRjlWrCEitHA2QYklTIhK5NaS0FCG0pGBLPr8/ZiadTGb2npk9l7VmXs/H4/v4fr+zb2vPzF7z3mutWWPuLgAAAPRrzdgFAAAAmANCFwAAwAAIXQAAAAMgdAEAAAyA0AUAADAAQhcAAMAACF1Aj8zs02Z2bk/PfY2Z/UQfzw0A6B6hC8gxs5PN7Coz22dmbmbrxy6TJJnZJWb2J2OXoy0z+3kz+7yZ7TezS83s6Ir7rU/X+x25n1/N3X50+vivpM/3wgaPvbFw291m9pe52+9lZv8r3fa3m9m/mNnxdd5D+trvN7MDZvav+aBtZpvM7N1mdpuZHTExYqFMd5jZN8zsd3K3n5M+54H0NR5cc31sNbP3mNl/mNkXzOzPzOzk3O2XmNldhdf+ltztbmZfzd32R4XX/QMzuyV9/r80s1PLtimAexC6gMMdlPS3kn547IJMhZl9r6QXSTpH0npJ3yLp5Usedry73zf9+Y3c8ksknSnpwZK+W9Ivm9l5dR7r7o/Mlku6n6TPSvqz3ONeLuk7JP0XSd8s6VmSvlbzPbxd0r9IOlHSSyW9w8xOSm+7S9KVkp5b9kZzZb2vpAdKujMrl5mtk/ROSb8q6f6Sdkj605rr4wRJ29LyPljS7ZIuK7z8n+Zf390/Wbj90bnb8q2qF6frabOkUyR9WdLvCMBi7s4PP6P8SPofku7I/Xxd0jVLHnOspN+S9BlJ+yV9QNKxPZTtKEkuaX3Dxz0rLdsXlRx8Py3p3PS2NUoO3P+e3n6lpPunt61PX+8iSfsk3SzpF9LbzpP0n0oO3ndI+mi6/BpJvyHp/yg5oP6dpHU1yvhfJf2TkgPlTZKe0/N2fpukV+b+P0fS5yvum62Hoypu/5yk78n9/xuSrqjz2MLzPCFdl/dJ/z8h/f+hTd+DpIel++79crf/o6SfLjzHhqTKXViuCyR9UpKl/18k6Z9yt99HSSj71mXro+S5HyPp9tz/l0j6kwVlcUkbKm57g6RX5/7/b5I+3ud+xA8/U/ihpQujcfdD37KVfFv+pJIWg0VeI+mxSlok7i/pl5W0Th3GzM4wsy8v+PnRjt+OzOwRSg5Gz1Lyfk6UdFruLi+Q9DQlB/xTJH1J0u8Vnua7lbRcfI+kF5nZue7+t5JeqXtaJR6du/+PSrpQ0gMkfZOkX1xSxjMk/Y2SVomTJJ0l6fqK+/7+gvW3c8nqyHukpI/m/v+opAea2YkLHvMZM9trZpelrT0ysxOUrLficz1y2WNLXCDpHe7+1fT/R0m6W9LT0266fzOz59V8D4+U9El3v31Jueq4QNJb3D3rhjzsddPy/rukRzZYH5nvknRjYdkPpN2DN5rZz5Q85tp0fbzTDu9qf5Okx5vZKWa2VtKPKdmvACxA6MLozGyNkpaEa9z9jUvu9+OSLnb3z7n7N9z9n9z968X7uvtn3f34BT9v6+GtPF3SX7n7tWmZflWHB8KfkvRSd9+b3n6JkoP8Ubn7vNzdv+ruNyjpCnrmkte8zN3/zd3vVNJydtaS+/+YpPe6+9vd/S53/6K7l4Yud/+fC9bf5iWvk3dfJa2Smezv+5Xc9zZJ366kO+yx6X3emnue/OOzv+9X47GHpCHh6ZIuzy0+TdJxSlqtHpLefomZPbnGeyjeVixXLWkgfoKkN+cWL3ruZesj/9ybJf2apF/KLb5S0kYl4fsnJf2ameX3tycoaT38ViWtr3+V21f/TUn37OckfSV9nl+v8TaBWSN0IQSvUHKgeMGS+62TdIySb/ohOkVJd52kQ60SX8zd/mBJ78paiyTtlvQNJeN4Mjfl/v5M+pyLfD739wHdcyCucrp6XH9m9mO5gddZy8cdSsZIZbK/b1eBu9/h7jvc/W53v0XS8yV9j5l9c/o8+cdnf99e47F5PyTpPyT9Q27ZnenvX3f3O919p6QrJD2lxnso3nZYuRp4tqQPuPuncssWPffC9ZExsw1KWqEudvd/zJa7+8fcfV/25UXS65WEzez2a939P939y0rGcD1ESbiSkhbdY5S05t5HybgzWrqAJQhdGJWZPUNJa87T3f2uJXe/TcnA5ofWeN4zSs4Ky//8WAfFL7pZSajJyrBWyUEpc5Ok7yu0GB3j7p/L3ef03N9nKGlhkJLxNV24STXWnySlZ6dVrb9iN1VSSPe3+j0Dr78vXXyjpHyX6KMl3eLuXzzyGY58yqw47v4lJeu4+FylZck/trC82IUnSVl3adV6XvQebpT0LWZ2v8LtVeWq8mwd3sp1xOua2X2UbL8b66yP9EzH90r6DXf/4yWv7zpyXVXd/mhJl7v7f6Sttr8j6ewF3bkAJAbS8zPej6Rvk/QFSWc1eMzvSXqfkhageyk5g+rojst1jJJv7y7p4ZKOyd12iSoG+ysZS3OHkoHq36Rk/Nndumcg/c8rGfz+4PT/kySdn/69Pn29t0pamz7XrUoHSUv6aSUnDazJvd41kn4i9/9zlLSULHpvZyhpCfnvSk4WOLHJ+m+5Ps9T0iL3CCUD1v9e0qsq7vu4dJ2vScv2p5Len7v9VUpaqE5Q0u11s6Tz6jw2vc9p6TY5YsC8pGslvVHS0UpadG6VdE6d9yBpe7q9j5H0g0pOUjgpvc3S5Y9It/ExxX1WyRjFryo3GD+3j+xXcjbtMZJ+U9L2muvjVCWtmr9Usa7PTx9nks5W0lV4QW5fPkvJZ+y+kv4/SR+XdO/09ssk/bmSLtl7S3qJpM/1uR/xw88UfkYvAD/z/VESYO7W4Wcw/s2SxxybHgA+lx6MrlXHZy+mB8bDfnK3vUnSKxY89gIlY12qzl58YXrwuj09IL4yvW29Dj978fOSfjn3vCcqCV1fkvSRdNk1ahi60vt9p6QPKhmLc1N2oO15W79Q0i3pa16WDx1KWmZ+LP37mZI+lQaQmyW9RdKDcvc9WtKl6fPcIumFudsWPja9z4sl/WNFGU9VMl3IHUpO6vipBu9hfbo97ky377mF24r71KcLz/1GSX9cUa5zJf1r+tzXKHdG7ZL18bL0tfKfrztyt7893U/vSJ//BbnbnpS+j68qCZ//W9KZhf3xreltX073zbP73o/44Sf2n+y0ZAA1mNn1Slo/6nSNNXne9UoCw73d/e4unxsAEIajlt8FQMbdl50dCABAKQbSIzh25OVa+hz8PjmFMwiXDn4HgDkws9MtuZTW7vQ4c3Hh9l9ML3+VzQ1oZvbbZrbHzHaa2WNWLgPdiwAAYOosufboye7+kfRs4+skPc3dP2Zmp0v6IyUnpDzW3W8zs6dI+lklU8c8TtLr3f1xq5SBli4AADB57n6zu38k/ft2JXMlZhdqf52SK5zkW6LOVzq9jLtvl3S85S4a30bUY7rWrVvn69evH7sYQbvzrm+MXYRZOvbe9zpy4de/evj/R99nmMJgGHW2b/E+Je5cc0xHBUKm9POII1x33XW3uftJy+/ZndPtWP/akVdya+U2/eeNSi9Sn9rm7tvK7puevPRtkj5oZk9VMuXJR80Om6ruVB0+YfXedNnNbcsYdehav369duzYMXYxgrdzX/EqIujb5lOOO2LZwT3bj1i2ZsPWIYqDARS3b3Hblm3/MrvWblx+J9RW9llEOTP7zNCv+TUd1A9rpcajQ96oz3zN3bcsu5+Z3VfJPHM/p2Taopcqud7tEXctWbbSmKyoQxcwtqxCrxNs12zYesSBN/8/AWxa6oasok0HdhO8gJ6Y2b2VBK63uvs7zexRSi5xlbVynSbpI2Z2tpKWrfxVQk7TPVcJaYUxXUBL+W/QxW/TtC7OS52ARagGxmVJqnqTpN3u/lpJcvcb3P0B7r7e3dcrCVqPcffPS7pK0rPTsxi3Strv7q27FiVaugCgtratV1ngKmvtLLPpwG5JdDUCHXu8pGdJuiGd6FqSXuLuV1fc/2olZy7ukXRA0oWrFoDQBfSgahxJ3YMuwrHq9iq2cDXZB+hqrJb/jC1qWS7exhiv+XL3D2jxRd2VtnZlf7uk53VZBkLXDDQZd4T6du7b36oCLzvo0vUUl7aD5POPX/aY7DU2l9w298/yKsGp7ecW6AKha0Y2n3Lc7CvrEBG4wlU3TOW3YRcnRyx7XN1WHpQjeGEsDKQHOpJV4ssqc7oX47Dq+K0+nrsM4aGdnfv2H/YDDIGWrpmhtese2WDlvFXHzyw6ANKlOA2rbrc+QvfcPtfFlqpV3z/BFUMhdM3Qql0TMXdtlAWtstvbhq9sfRQrcebjmoYhAlfb14j5c9lG18ELGALdizO3+ZTjDv3Uvf+i/0O16cDupYGrzf3zlXzV3wSu+MR+9YBYPpddm+v7Rjxo6cIhU/ym2CRoVT1+WatX1TrbuW//Ya8f00F7quY0nm6Kn+c6mr5vghqGREsXDtO0Agq5Ul81cK36PF29ProzVPDNXqdsHN+iHzRTVf8QpBAqQheOULfCCjFwZd2CXQeeVZ+PA2r8hmglYz9pjuCFmNC9iEkYolWp7SB7DqThq5pnq+o+sZhLF2PVvFtly4rrgzm7MCRauhCtvlq16rxuXTEeqKeu7LI8i/6vWoaw1A2XBCyMiZYulKr6hhjCt8IQxkr1dU28NjOgo7ll64/1G6cQ6idgEUIXGhujuyKEoFXExYjnbdG0EkOdJVk1L1zx9rxloST27kiCF0JG9yKCF2Lgyiwr25ymKMA9hmgpqxOO8vPw1Z2PbwqBZdG6iT1UIm60dCFoIQeuzLIWr64nR6Xra3x1Jk/tejtxll4zxRYvwhZCQOhCsGIIXJm6XY0H92xf6WBM4BpfmxC9bLtXtYiO1X3d9VmPY4Wfuq2BwFDoXkRtVE6L1T2bsk2XI5NnhmGV7uKmj10WuJpcvmtMsV46DOgDLV1oZC7z/qyqy+BF2Brf0NdiLAtcQ4eVLlqnCFjA4QhdQMAIXGFqul2adC0WA1cIwYUvW0A3CF1oZMiKd9fajVGN6wKWWdTCGULrVpdiLjvQF8Z0oTa+6Q6HMVzTVLVNme8NmAdaulAqlIA1x9Yuwta0rdmw9bAWrxgCVyj1ARA7WroQvBgOSl0hcCE0bQMXQQ04EqELpUIbj7Fr7cbJhy8C1/xMfZ/OrtdaXAbMFaELlUILXtL0D1II36rhuMl0IVloiT2ohPw+QiwTpovQhYVCDV5TCl8Mmo/HkIGriHDQjxDrOEwXA+kRrSx4xTzQnrAVB7YTgC7Q0oWlQv8mOKVWLwDAdNHShVqKwWtZV8fQF7id49QSCMOqFzFH/5rWX0BfaOlCK01av7IL8/bdYhZbixcH6vhlY7TqjtVimw8jX+eU1Tuht95jumjpAkbAwTd+xaBVt8WrODkqmstCU77FqmmQInhhDLR0AQMjcMWvq9BEl3hz+bCU/U2AQiwIXejFWPPyxNbFiPiM1UpFsChfB6wXxITuRfRujOAVWgsCrVtYBcGi+3WQ1UusWwyJ0IVWQj/7J9/iNVYAI2jNA9u5f30FruxvgheGQujC5LXpclw1qHEgnqZi12IX23nTgd10iwMzQehCK2VnDxVvy4TeKlamSUsZAWt6OLswEcJnt49WKFq2MBZCF1ZSnAS1aqBrCJV3W4SquFUFqLG2a5tAN0YXWAifWcIRpoazF9GZRRXkUBOkdmXTgd3adGB3NOVFub5brJo+Pxe8BuaNli4MLvSWr9DOfER9Y3QL5l+TVlEAixC6gFQxbHEAjUuXgavNrPF97y9DtrqG/KUIiBmhCxCBK3ZjDnxvsq80DXN0bwPTQujCKEI+w5HAFZe2gatpt2D+Ptlj2+wrXHsRmC8G0iMIMQ2yRzi6Ci/LnqcYrtZs2DrpcM5nEegHoQtBobLHWIZsfZpyYANQrbfQZWanm9n7zWy3md1oZhenyy8xs8+Z2fXpz1Nyj3mxme0xs4+b2ff2VTaEbejgxdmKyBzcs/3QDwB0rc8xXXdL+gV3/4iZ3U/SdWb2nvS217n7a/J3NrNHSHqGpEdKOkXSe83sYe7+jR7LiEAVJ11ddPui+9UR4gWysdyQwWiVMVxtjdnqG9IYS2BKegtd7n6zpJvTv283s92STl3wkPMlXeHuX5f0KTPbI+lsSf/cVxkRh7oHn0WXJgK6cHDP9kGCF93swDQNMqbLzNZL+jZJH0wXPd/MdprZpWZ2QrrsVEk35R62VyUhzcwuMrMdZrbjC1/4Qo+lRqw4YKFPXbSwhdx9yZcWoD+9hy4zu6+kP5f0c+7+FUlvkPRQSWcpaQn7reyuJQ/3Ixa4b3P3Le6+5aSTTuqp1Ihd/rJDdUJY/gLXwDKrhCYCFzBfvc7TZWb3VhK43uru75Qkd78ld/sfSvqr9N+9kk7PPfw0Sfv6LB/mo8mlh4bqQkLYivvAVIOWRNgChtLn2Ysm6U2Sdrv7a3PLT87d7Qcl7Ur/vkrSM8zsaDN7iKQzJX2or/JhfpgLDHWVhe6+5uYqtrIOvY8SuIDh9NnS9XhJz5J0g5ldny57iaRnmtlZSroOPy3ppyTJ3W80syslfUzJmY/P48xF9GFZqxetXHEYa2b37HVjbBElYAHj6vPsxQ+ofJzW1Qse8wpJr+irTECm6kzH2A6iGFed4LUoGJaNJdy5bz8tsmrW4pd9jllvCB0z0mPW8pU0gQttVIWqRZOs7lq7MdiTNxZ1ww8VatoEruLfQIgIXYCSWelDH+yMcHTdujWEOoEkH3aKZwBnt9GahFil01Tdama7Cst/Nr0Szo1m9urc8s6vkkPoAnIIXlhVSIFr5779h366FFLwKpaF1i4scLmk8/ILzOy7lUzOvtndHynpNeny/FVyzpP0+2Z2r1UL0OuUEUAImo73qDpo0v04D3W3c34gf52wPlTgWiV0hDCerIvXD+F9IDzufm06WXvez0h6VXo1HLn7renyXq6SQ+jCpJWN9yhWxsn/i8+EI3CFKdsuXbVQNt3Odc+gjCFw5Z8jtsBSdkZy/v/Y3s+cPPDEY3XxD2zu5LneePln1pnZjtyibe6+bcnDHibpO83sFZK+JukX3f3DSq6Ik/9wl14lpym6FzFZVQegquVVB1wCV/i62EZ9bedQB8xPDcEKkm7LrliT/iwLXFLS+HSCpK2SfknSlek8o7WuktMUoQuzkw2aL55dVtZiQeCKR4jbqm3gatNiNeRYpr5ea9XgRPBCC3slvdMTH5J0UNI69XSVHLoXMSubDuw+YtkYY7jqdoeFGCRCNtaEqX1o0s03VOBikDom6H9LepKka8zsYZK+SdJtSq6S8zYze62kU9TRVXIIXZisJtdbHNKiuZsyWTjM35cAVs+Uglcdfezjqz5n1eTDY4lxnBq6Z2Zvl/RESevMbK+kl0m6VNKl6TQS/ynpAnd3Sb1cJYfQhVnZtXZjaWvXUIphoKr7KVueL2uMl52JSRfrtziwf9OB3bMb0xViuAmxTBieuz+z4qb/q+L+nV8lhzFdmKxVv2X32VpSd0by4v0WzXKOceUDW1fhOJSWoiZWnbKi69cmcCEktHRhVF1eM61uhT1ka1dZQGrT8lHseqTbcbEhuxj7Xv8xdo3FGBaBIdDShdEU59BqW1G3eWzd4BNiq1Kx7CGWMQRTCqNV+zjhBogLoQuj6OpgMcRBp8tQ09f4HoJXuSbBa0ohLWarfKZjaxHE/BC60Ngq13Nb9LimFeaqgatJAAot1JSVPbQyhmLNhq2jBKrsNbvsyi67wsJctK135raeEDbGdKGRRZfaKGoSopbdt6+Ks8n4rizUhNIiMvaZmLFZdMmgULZpHXMIETGOYwPqIHShtr6+YZZVrkMeWJqGl7qTqY7R8hRaMIwB66p/+c943c928dqJhDBMAd2LWGqVQe7LhFKRrjrWqk7g6mM8V9Vz0tVYD4Grf+UXmO/HHFoBETdaurBQn5VYKIErUzYhaR1lB+65TYgZoyED1xwnSV2kj6tFFFvGgBDR0oVKcwpceU0OjmUH7qFbmTiYY+4IXIgFoQtH6LM7UZp2pThWtx5nM4YpH8o56aG9RfURgQsxIXRhULFUim1aj0IMOSGWaW4IXt1YNNUMA+0RC0IXjtBX5RVbpVgneGWhJoRww6D6cM05eHU1Lx8wBYQulBqzQgzpDKQmwavJY4ZG8BofwQsAoQuVxmi2D7FybhKixg5cY79+DJgmYhj5+qOqDulyAmUgBoQu1MKYiWq71m489BOCUMoRsrGC11wCX9WEx2VfqurUK9Q9mApCFxphYsN7hBS0ikItF6YfvNrUEYseQ+DClBC6EISYAlfIYSsvhjICGcIV5oAZ6TG6GAJXrAEmf13Jg3u2T76VBdNCEMPU0NKFUcUQuGKXD4ycxYiu9XUmJoELU0TowmgIXOMgeKErWeCqE7yanIxD4MJUEbqAGSh2jxK80JeyL1NNLi1G4MKUEbqAmSB4hWXKE6RmIavpdVwJXJg6BtJjMHQnji8/sF5icH1XFgXYKa3f/L4Ty1m8Q9m5bz+hEUsRugBgYGs2bI2mpbGqRa64vE4Am2ow4Qsl6iJ0oREql260OWB1pdjahf4sauVas2GrNu3Z3uu2L9vOxdcbel+YavAC6iB0oVKxciRwtbfswLbs9rID86LH0O0Tj00HdrfaXm3DUl8hq+37iEW+/iM0oi1CF0oRsFbX5cGt6XMVx94UMWlqWLJtUSe0hNxKWed9EFgwZ5y9CHRo04Hdh35CUVamkMo3ZU3HbZVtl/y2i2W7xVLOtvJfSvmCiiYIXThCVolMveLsSohBq0os5ZySquBVtTwfsGIOylXljTGk1J17LMb3hmHRvYjDFAPXXLqeYjugIS6xnKk4lKyemWJX45TfG1ZH6MIhVd/SYv72RpjCELIvJoSrwy0bXM+ZjJgbQhckVQerNRu2SpGELgIWpmJK+/JUz2pcNu0LLV4oQ+hCpezb++ZTjgu2tWtKBycA4alb9001XKJbhC4coSzIDB28CFPDmcOYPYxnjCAyVisTwQvLELpQO0xlFVhf4YughVAVx2oVgypjucr1FUDqTlTa15ix5DkP3wdC7Q1AWAhdM1VVQdQJPl22ehG0ECNC1nJZ4Fo19LSpa/J1VNfBq+q5GLuFOghdM1SsxMqCz6Iup1UCFyErPHOZFqQtAlZzXQWuRZY9d1evTZhClwhdM5EdOFa52O2hg0+DLgNCFoCulV0Tdlk4YnoKhIDQNQP5b+rZQM9WYavwHEUELEwNrVzNjTWQnFCFGBC6ZmhZOKrT1UTAwpQRtsLQdCjDolavbJxXVTCrOzgfWAXXXpw4Dh5AM3xmwnBwz/bWX+4WhTXOMsSYCF0RyS6wWrfS4OABIDYH92w/YkhEG1yMGiEidEWirAJpEsAALMcXlTAVg1fdeq9tPUm9ir4QuiK26cBuxlahE4QNhKxscH6dMLX5lONKx2cRqjAWQlckqioPqfqA2WbuJWbaxlyxr0/HovoSGFNvocvMTjez95vZbjO70cwuTpff38zeY2afSH+fkC43M/ttM9tjZjvN7DF9lS1mm085rrSFq4sDBoFr3tjeWFUX00XsWrvxsOdp+px1w9ay1i5aw9CHPqeMuFvSL7j7R8zsfpKuM7P3SHqOpPe5+6vM7EWSXiTpVyR9n6Qz05/HSXpD+hs1tZ1ZnNnIkZnr7PQEzjDkg05Z2Koz3UNXt2W302KGLvUWutz9Zkk3p3/fbma7JZ0q6XxJT0zv9mZJ1ygJXedLeou7u6TtZna8mZ2cPg96MscDLIBhNJ37qouA03ULVd3nI5yhjkEmRzWz9ZK+TdIHJT0wC1LufrOZPSC926mSbso9bG+67LDQZWYXSbpIks4444xeyx2iVb+RE7KwzFxbu9CvPluN8he4rlIckjHWzPmYt94H0pvZfSX9uaSfc/evLLpryTI/YoH7Nnff4u5bTjrppK6KOQscSLHIXA9CdC1Ow6JAx1neCEWvLV1mdm8lgeut7v7OdPEtWbehmZ0s6dZ0+V5Jp+cefpqkfX2WL0ZNr5soNQtbHIDmK39Nzc0jl2Uo7O8AhtTn2Ysm6U2Sdrv7a3M3XSXpgvTvCyT9RW75s9OzGLdK2s94riMta57PZnPODiZ1A1dxFmjM2859+ye/T0z5vY1l0RfCKY95mvJ7Q7f67F58vKRnSXqSmV2f/jxF0qskPdnMPiHpyen/knS1pE9K2iPpDyX9zx7LNgt1DyocfLDIFE+dZ58fTshzZnVRtlDfG8LU59mLH1D5OC1JOqfk/i7peX2VZ0radDHmccBBlXwXY36M15QG17P/9+uwbuqBA0mdAfVSsm/ny5b9PcUvGDicmf28pJ9QMmb8BkkXSjpZ0hWS7i/pI5Ke5e7/2cfrD3L2Irq1+ZTjdHBP88dxsMHc8RmYvmLQ27lv/xEniVSFQcLXtJnZqZJeIOkR7n6nmV0p6RmSniLpde5+hZn9gaTnKpkrtHNcBihSTVsdONigLs70wirGauWq0qYLseox2fL8D6JzlKRjzewoSWuVTEv1JEnvSG9/s6Sn9fnimDDCFtooC15T6mLE/LQJSISq6Kwzsx25/7e5+7bsH3f/nJm9RtJnJd0p6e8kXSfpy+5+d3q3bI7QXhC6AADAKI45+WQ9/KUv6ebJLv/L29x9S9XN6bWez5f0EElflvRnSi5BWHTEHKFdoXsRQG0htpwWp0kBgArnSvqUu3/B3e+S9E5J3yHp+LS7Uep5jlBauiLFAQZjGbqbscm+nr8vXaHjoVsOgfqspK1mtlZJ9+I5knZIer+kpys5gzE/f2jnaOmKEIELYxtqH2Rfj8uaDVsJXCvizMn+uPsHlQyY/4iS6SLWSNom6VckvdDM9kg6UcnE7r2gpQtAkPoIXGs2bCXI9YSWxW4QWvvl7i+T9LLC4k9KOnt7w6ZkAAAgAElEQVSI16elC0AroYeX0MsHLEKL1zQRuiLDgQQhCX1/DL18U0ErVz927tt/6AfTQOgCsJLQg03o5YsdgWsYBK9pIHRFZs2GrVRymLyug1L2fASwblEX9ac4touxXtPAQPpIMSAYIVtlWom+9ms+L+hT1hLVZTjKLuC96DmbtIAR3MZHS1cEOFggdPl9lFaleaCVq1zX3YDLglKTIEUX5fho6QocBy7EomxfHfN6jYtel88Vupa1Skla2jrV52svM3TZcDhaugK26MDAQQOotizoMTayvV1rN2rX2o1jFyNIY4YZWrziQEtXIHbu269NB3bXui+BCzEZo7Wr7msyNrIeQlZ9TVqdukSQigOhS/fsrPnQM9RBouy1y4zZTQOsqsn+WzcIFZ+v+Jiq5yg+Lvuf8HUkwlYc2gQuuhnHQfdiTr7ZvI8KeNVvIgf3bOfAECHCcnN1ugebPiaz7HOUdT3OtQsyqwcJXKsZquVplWMCk68Ob/YtXWU726HKpsNvAst26joVXNYaRpdIPLKDNtusuap1VhWEmq7ffOtb9tiqMDf1bUfA6sdQrUmbDuxmG0ZitqGrLASVdfEd3BNOS8WutRsPfYDncCCYguzAzrZq10Xe97qb63bhAN2vocZ15fffRUNU2N7hmGXoqhu4+nrtLDjRpDsPcz2wd2XILz0xjJ3kABqHLHj11drVpF7Jjm/sO+ObXegaM3AtKkMdDHpE7PIHilDDTYghmYNl3LoMXqvsn3RDjm9WoWvMwLVqc3PZB3bX2o2tv8EMHTSXlW/o8mB8MbQqjYUD4/R0Eby6+EJAq9e4ZhO6ioFnjIN8H/38bT84oX3gCGXzFEPL16qq9u3iPh3aZ3IRWt2H1Ufra7HVi+Euw5hF6Ko7F1aVJmdQLTPWxHmxC/WA1PTASXistmwah2W6OquxS4v2h1D3aXSni/q+z/2X7sbhTT50NZnpfShtPohMZBemphVWmwoutP13DIumdAgZBzW0qbdDHFeIbkw+dE3pgEXwmqdFB+0p7d91MA4MU0fgmrbJh66pyQcvuipRFsimHsSaBC8OYIgF++o8ELqAiak7cDtmxQMUrV+IUYjDX9CvqEPXnXd9Y+l98tMq9KnvMSdV3Yq0dqGuObaKjW2q47kY6rC6RfU2V7GYrqhDlzSPMR5UbujLVFrFmtQD+fecf591Ht/0QDjlgfQEr9Xcs+7K97uhgteU99EQRR+6lhmylSv7e1nl3aRlikoNY6hbCYcUzhZ99qreT9YSXjewLbofLROIVUif46mLPnTtWrtRm8cuREFZ5d9nFyBdjBhLX9+Qlx0Eql63WBdkrQWLvs131VI+xy4hWruAZqIOXcce/NrS+4xVEfbR7VkMVvnKjuCFKWkb5so+d/nglb/tnnqhu89pWX0z9e4bghdQ35qxC7CqTQd21w5VfVV8Qw+ezywKYcBcVV1BoiyM9fHZLXvOTQd2l/5MBV/4gHqiD111ZYFryG+cWeXfRYVUFrAIWUC5sbv56oa5KQUvdGvsfRj9mEzoyu+gy0LOkMGrr8AF4HDFADP2QYvgBaAo6jFdRVUBZ82GrVLhtlXn76ozZmvK4ziAkFR9lquC19SnmQEQpuhbunat3Xjop81jV3ndPtGaBdQTckvRVMNd2TpnXBewXNQtXXeuOWbl5+hqlu41G7Z2VuksClyEMeAeTT+rQ12hIq/OGdQxneGYrb/sd77cnMnYjSG7xmPZ76Yi+pauPixrOaucI4jKBhhMk/BU/EyPPd6rTMgtdhjOUPtm2x4irCbqlq5F8t+42rZAtflW3Od8WcVQV/xWSfN+eBbtP1R47bRp3Rpb3fkCQ2/xKlv3xTLT2jWOkPcb3GOyoUvqJoQsCl5DVS6LXoOgNZyuWyLaPF9WsdZ57JQq4bbrPqR1UBzfVRXCyrrtQtdH8MrXbYS4cjHtI0hMOnR1JduxD33wlwSdIWaHJ2x1L4bunSZljD2Yrbo9Qn5v0vLWr9BbvYq6LC/122Ix7Rc4HKGrppAuuUOF1E4MoWpoZetk7Aq9i+009nuYq3zwatvatXPf/ihb+4bCOokboaulLHh13cVY9VwEreYIWe3k19vQFfzcAtfUWrvqyuqzsnGq+X0ge/9zGSe2bNzfFPeFueHsxZoWhZ6+p4ogcDUztevajWnI9Ti3wJWJZS6vuttn0f2yL6rF/w/u2X5E4Cp7LBA7WroALFTW2pJf1kXL2FwD19RVnUVedRZkk+cDYkTowqTQwtW/4uSYZbdJw4aguq8V6v5Rd0qJsayy3paFrVWeb0oBLOTtj+7QvdhAkwtPN60MplR5YFr6br1a5UA8hwkeQw2KixTLTJc/kCB0rWjzKcetHJi6mIeL0IY+LLsyw6otTKsGLoRrUYtoV+Y2zovwGj+6FycgC1xlU1lUdfeEOFUAhlF1vdG62794v7LxXsXlXe1b7KPDW+UgT0Dox1TPap0DQteElTXx170/H+jpqLMtu9zefe07XT5vLGcMzh31UDXmMosT3YsNdd2cvaxbcKxuw6Yzn9PsPb6su6/4E7K2rWttxLx/jln20PehKVhlEH1X+0bMn4+Y0NLVQlenLYc+DqvON6lFrWlU1v2ay/qdy/vEPHVx1mKxrq4zfKT4GD5nwyB0TUBV69uii3U3URW+mnRXluFDXk9+O85tnXX1fqfwLb7qQDrE+J6u6hL0a9E2ytchZbP+s32H0VvoMrNLJX2/pFvdfVO67BJJPynpC+ndXuLuV6e3vVjScyV9Q9IL3P3dy17j2INf66Hk/Qu9hatK1YSYXTxfV2KrPKpOblg2WH1Khn5vUxrPVTxDsM91OfRnq84JGmiu7cS06EafLV2XS/pdSW8pLH+du78mv8DMHiHpGZIeKekUSe81s4e5+zeWvUhxZwn9Axlr4Mp0/eGc01k4y97nXNbDkGI6mHQxQWrfn6eqM1/7VmyZOVy8IZoJUeent9Dl7tea2fqadz9f0hXu/nVJnzKzPZLOlvTPTV93qDFFQ16OomwqiLpia/mZsq63w1xCWtvP9Fz3+6G/yPRVx9DthSkaY0zX883s2ZJ2SPoFd/+SpFMl5SP/3nTZEczsIkkXSdIZpzxo4QuN2Rzd91mOxUthxDxJYJcHiTlV0HXfa9W6rXp8SGGu7ESNOuWb035QZkrBa6po5ZqnytBlZmsk7czGY3XkDZJ+Q5Knv39L0o9LspL7etkTuPs2SdskacujNpbep2iKZ9QVQ1gxeOX/n8u3xTm8xzaarpe267HOftZFS1UxULDdy9VZL1OpD4FYVIYudz9oZh81szPc/bNdvJi735L9bWZ/KOmv0n/3Sjo9d9fTJO3r4jWLploRlV0ENqbgFXLZUE+dbbjoVPamc8O1+ZxOaRB9F1btDeBz2w6tXOMxs/MkvV7SvST9kbu/asjXX9a9eLKkG83sQ5K+mi1096e2eTEzO9ndb07//UFJu9K/r5L0NjN7rZKB9GdK+lCb1+hCF92SXXb5LRs/Vha4yoQevDBPbffJpo8jcFWLtV7YtXajNo9dCETDzO4l6fckPVlJY8+Hzewqd//YUGVYFrpe3vaJzeztkp4oaZ2Z7ZX0MklPNLOzlHQdflrST0mSu99oZldK+pikuyU9r86Zi30L4RTlsuC2SpgjeCFmWXCipQBS/GeDY3BnS9rj7p+UJDO7QsmJfGGELnf/h0W3m9k/u/t/qXjsM0sWv2nBa71C0isWvd5YqroyhjyDsa2yFjeCF2LXxfQKAGbnVEk35f7fK+lxQxZg1bMXj+mkFBFoOoYkpLMJCV5AHAiTmJs71xzTZW/SOjPbkft/W3ryXab2SXt9WTV0DVrYscU8kWe+RS4/uF6KdzwHAMSGUN2r29x9y4LbBztpr8qaIV8M3dh8ynErdWsWHxtrkASmaM2GrQz6nygC1+g+LOlMM3uImX2TkivhXDVkARaGLjN7vpmdsOguHZcneHVbhapCUZfdjl2OJ9u1diPhC5MVY4iJscxAyNz9bknPl/RuSbslXenuNw5ZhmUtXQ9SckrllWZ2npkVQ9azeipXNIYYuxX6YH1gaHNpMSB4Ad1y96vd/WHu/tD0BL5BLQxd7v5/K5kz602SniPpE2b2SjN7aHr7rgUPn6xia9eYg+azrsZVuxwztHYhZAf3bJ9N4MoQvIDpWDqmy91d0ufTn7slnSDpHWb26p7LFrRNB3Yv7WpcNQgNFeaKoY3gBYSF4BW/uX1ZQLllY7peYGbXSXq1pP8j6VHu/jOSHivphwcoX/DqjPGKqXswprICmLeQpuYB6lg2ZcQ6ST/k7p/JL0yvy/j9/RULdRQnZ11WAdUNVMm8XszjBQBAl5aN6fq1YuDK3cYROTVma9fOffsP/bRRVa7NpxxHNyMQELoY40XXIjKrTo6KFoa4fFD2/FVhrM7r0+IFhKUseHFAB+LB5KgdWVbxFcNPH2MRys5iXDXc0eIFhG3OLWB8IURsCF0jWqVbsGhRuOpiBnuCF4AQMZgeMSF0BaCL8JV/fNk4r1Vfg+AFhGvOrV1ATAhdAekieJV1YxbDV1tMJwEgNHQxIiaErg51MaA1C0V1wlHbELRK8KrzjZpv3cDw+NwB4SN0BWiIMQp9vUZW8XMAADCUkM/gDLlsGB6ha8baBi8CFRAmPpthIXChiNDVsdg+ZF22eOUr/NjWAzAVBC8gXISuHsQWOFYNXms2bKWiBwLC5xEIE6GrJ10Gr/zg+tDmpKFyB8LEZxMID6GrR10GpLJpH8YMYIsq9Nha+gDEj3oHMSB09Sg/f8yqs8IXhdLiVazoyio+uh8BDCGUehGowgWvB1YWvEKoKFa5CPeib5iELYTg4J7t7IsARkdLVwCGmOl9rNnks0BG0z8AYO4IXT2rGzb6DEXZc48dvIAxsR8CGBuhawB1uw/7DkWLytHmtemuQWwIXgDGxJiunu1au1HS4YFnUcDpY7B9V2PGysZ9rdmw9dDzc+FZAACq0dLVoyxwFQ01cL5JgFtWpjplrnq/ADBH9AagiJaujtUNHmUhpuvuxabhLrv/onKUtXZtPuW4IM7ABDA9WZ1apyU9uU9YQWfNhq10a+MQQtcKum7ZqTvmqqqrsstuxOJzF28neCFWc5k+ggM9EB5CV01jd50tm0dr6MBTVZ5dazdGMbZr0faMofzAIlMJXGPXu0DXCF0VQvywlwWrMVuX8sErplauZds2luAIlJlK4Mpr8nlcZaJnoG+EroIQw1ZIipXfwT2HrzPCCmIwly7G2FGfYGoIXamph60+Lz8US8XYZBvT2gXEKcTB9EBm1qFr6kFLWj4nWEzdgm213c5Nzppq+zoEO3Rpil2LTYXYglm3bsjqa7bjdM0idM0hXLU1teDVx7buc/9Z5bkJbKuZWhcjB+q45b8gM83EdE0mdBGs+pWFsxC63djWiSbrYextFqqpBS+EadmXWwb+z0fUoevONcdwAFb1eK2uP8ghBC+2dztMkTFtxVaRpl3jU0GAbo46dVhRhy5UaxK4mnQxjtkdSeXQj6r1OpcD9tRau/Lbc05nFse4Dbv6YkzdGA9CFyQ1D1P51q6qir2sIlh2O8KxbPtNyVSC16LP1Jy2Z2iKdWudsFXsUaC+nAZCV8T66j5c9tz5bsaiOhOPIl7Ltl/MB/EpBK9NB3Y3nhpl0XNhdW0CV4b6cnoIXRHqc9DlFAZ0Ft/DlM7ODF3sXZXZ2KjYwlcfZ7vFvi1D0Lbuoc6aLkIXgN7F1qISW6vXkNMLFLdliNtvTKsGJgLXtBG6MDlUWnEJdaxRLMErC1xjdUWFuv3GQN2DZQhdEQixyy9fJioadI2urbiNsf3GDsnUg6iD0BW4EANX0dRmtUe45jQFQl0xDbZu0irWZj7AsYMXsAyhK1AxhC1gTHOdADTv4J7tUkShq8yyaS7mvH0xPYSugBC0gOYIX8iLubWLHoPpWzN2ATANBEaMLaZuti7M5YLIbbZrqOGFehKELnSGCgVjm0vwGvuMxaHtWrtxEhd4DzUMYjh0L6JTDKrH2KY4DuiIkDWTsFVEVzJiR0sXOkeLF8Y2lRagg3u2z65Vq46mLV8x4MvqPBC60AuCF8YW+4E5P2Yr5vfRp1C2cZPAVHZfAtd8ELoATFrbA/OYZ8ARuJoZcx21CUz5xxC45oUxXQBmIZaxXgSudorravNI5aiLsDVPtHQBmI2YQkxMZQVQT2+hy8wuNbNbzWxXbtn9zew9ZvaJ9PcJ6XIzs982sz1mttPMHtNXuQDMG2FmHmhJQoj6bOm6XNJ5hWUvkvQ+dz9T0vvS/yXp+ySdmf5cJOkNPZYLADADfQcvgh2a6i10ufu1kv6jsPh8SW9O/36zpKfllr/FE9slHW9mJ/dVtlDxAQaGEXJrVyhn5E1FX/Uq9TXaGHpM1wPd/WZJSn8/IF1+qqSbcvfbmy47gpldZGY7zGzHl774xV4LC2C6FgWbsc5cPLhnexSD/eeOwDU9ZnaWmW03s+vTjHF2urzT4U+hDKS3kmVedkd33+buW9x9ywknnthzsQBMWUgtSkyCCozq1ZJe7u5nSfq19H+p4+FPQ4euW7Juw/T3renyvZJOz93vNEn7Bi5bEPgGBcwPgatfXdar1NGT5ZK+Of37ON2TQTod/jT0PF1XSbpA0qvS33+RW/58M7tC0uMk7c+6IQGgT7HM34XxEbi6d9OX7tQL33lDV0+3zsx25P7f5u7baj725yS928xeo6RB6jvS5VXDn1pllN5Cl5m9XdITlayEvZJepiRsXWlmz5X0WUk/kt79aklPkbRH0gFJF/ZVrjlYVjEMcYkeKiegnvxkqAgXdVoUbnP3LVU3mtl7JT2o5KaXSjpH0s+7+5+b2X+X9CZJ56rB8Kc6egtd7v7MipvOKbmvS3peX2WZkzoVw859+7k2IhAAZp+PQ1m9mrXOvPaHHjV0cdCSu59bdZuZvUXSxem/fybpj9K/Ox3+xGWAIhfit68QywRg3tp+2SRwzcY+SU+QdI2kJ0n6RLq80+FPhK5IrRpsaO0CwkEr1zCa1HtVdSyBa7J+UtLrzewoSV9Tcqai1PHwJ0JXRLpuQSo+XxchjFYuoB7OWBxH2zoqP9ibwDU97v4BSY8tWd7p8CdCV4BiDS6xlhvI62tiVAbMx4vAha6EMjkqAkBowtwNPRM9rVzhI3ChS7R04TCM9QKGQeCKA0ELXSJ0YTLyYZFWOwBAaOhexBGaBpad+/aPHnJonUPIiuO5aOUC5onQhVJ1Q9TYYQtYVd+XACJwAcjQvYhKBCpgNZyxCCCP0AVgtrpu5SJkAViE7kUA0DDTRdC1CMwbLV1opDhgfcwuyKrB83SLYgx1Wrk2HdhN8AJmjNCF2spCzuZTjhsl5JSVhbCFJrrqWpxbl2LVeiNMAssRuhAdpodAzLLQUhZS+j6Tsk+L3tdc0bKJIkIXokKXIroQQrgJoQx9mFvQWLYdu9jO+fVJS2PcCF2orewSQUOGHQIXEId8MFgWBprct0sxhd46ZW3yfgho4yF0oZGxAg5diuhKTAfbKWiyvttsm6l10w6hbP0QxIYRdeg69uDXSpcXdyh2prgtCly0cmFMcxtEHyICVjdYj8OIOnRJ7Zpd24SwOq9DuBsWgQtNdXlgIXABaCr60NXGsjEEbSvmZY8jlDVHtyJCROAC0MYsQ1fekE2qnHXSnUWtXJy6jqaGmI0eGMPHX/HK0uUPf+lLBi4JJEJXELoaPDoXxcBVtf4Y2wcML3+Q58A+rqrAVbyN7TQcQlekFgW1KYWLZeO2Vj0zakrrCsNZs2ErXYwFZQf4smUc4Pu3KGx1cX+0R+iaoLl0Y3bRNdz1xIWh6LrbPMT3OHfLDpRDhpsmB+3svoSv5ghH8SN0zchcWseGNsakhEOf3k1X7fj6bL0YIwARvuohaE0LoQuSxpsVuq1Y55Rp2sUZ6vuMbX+Rwl2X0vgH1jatZl2VuU4X5JwGg4+9L6Bf5u5jl6G1TeuO9z/7ge+a5AcvFKEdUEM+cM5daPuKVH9/aXr2Yh/juTjYrm7MFruiJmUZettvvOwqSdLuC5962PJHXP6X17n7liHL8s1nfKuf/ct/1Mlzve9nv3Pw8jc1iZYuzsLoT0gtGgSusM1l+/Q1gP7hL30JwWtFIa2/kMpSZeNlVx0RvNCvSYSuvOKOTgjrDvNfYaoO7tnOXF2Yjd0XPpXANZLJha4iWsG6tyx8zeXsScxP39NE0NqFoRC4xjH50JUX+tk8sWnanbRqV+WutRtn04WF8Aw1LxfBC5iuWYWuJtpUegS1+tp2VRbvTwhDiLKuSiZQBZBH6OpQ22+ncw5rjBNDKLoa19XFc9DaBUwToSsADP5PwhfBC1PQRXgjcAHTROgK0FwH/xO8ELIm11tcpVuRwAVMF6ErcLSCAeHgQtcAVkHoikydS2bE/HrLWrsYOA/M2+sv37nw9oufs3mgkgDNEbomYFF3RJuAtKx7Y+ggJq0Wtj7+ilf2sh6aPmfX2wndqzMeK+TWrqrLu0zBsrBVdj8CGEIziWsvIkxjBompjIshjI2rKoD1FbpW3W+Xha5YW4nqBq4qob6v0HDtxf7R0oXejHVCwFQCl3TPeyF8TV8X++2iFq46wSXEVqJVA1fZcxTf2+sv3xnM+8W0EbowiKHCw5QCVx7haxxl3Y1dtnKFvL+OEcC6CFhtX6dsGUEMXSN0YVBtx1fVed45KIavqvdNOOtOqF2JQ1rWUtTlc4ckxJY/xG12Y7r6rDxQXxehIKaD1lgIX+EYc3/tM9i0rUNDDlvLTPW4wZiu/s2mpavqA04IG0eb8V6ErOb6allEIvR9cohg06Y1KObAJVWXn+MHlpl86Gr64aZff3iLAljoB7UYMB6sH6HtmyEEmToD1qeM4weWmWzo6vLDzQdpOKEdyKZkrpeX6kNo+2moYSbUcg2J3hTkTTZ09Y0PEmLWZWiYW4ALKXARauKTbTOOGfNE6OoIIQxzNcYVCsYyduAiZE0HZ0bO02RD18XP2TxqBcW3GTQxtdA+tUseEbbQJwLYfEw2dEn37LyEL4RsbgfUJgEmhIBG4MKQCGDTNunQlSF8IUSL9kcuS5IYO/CMibAFAtj0zCJ0ZcbucpQIX0iMvR8iTOwXqMKxYxpmFbpCwgcIQB6BC3XQ+hW3NWMXYO5ef/lOKtuZYXsD6ALHj/jQ0hUIvr3MAxUkgK7RcxIPQleApjZ9ABIELgB9InyFb5TuRTP7tJndYGbXm9mOdNn9zew9ZvaJ9PcJY5QtRFkTMgfteLXZdmzv+WBbo0scL5ozsx8xsxvN7KCZbSncttnM/jm9/QYzOyZd/tj0/z1m9ttmZsteZ8wxXd/t7me5e/bmXiTpfe5+pqT3pf+jgA9TXNheWIb9A31h32pkl6QfknRtfqGZHSXpTyT9tLs/UtITJd2V3vwGSRdJOjP9OW/Zi4TUvXi+kjcjSW+WdI2kXxmrMKGLeQxYVUUQ2/so03Ulx3xdAFZBl2M97r5bkkoaq75H0k53/2h6vy+m9ztZ0je7+z+n/79F0tMk/c2i1xkrdLmkvzMzl/RGd98m6YHufrMkufvNZvaAsgea2UVKkqVOvs+xQ5U3aFM5MMdQOYzxzXEq23fONl521aG/d1/41BFLgrmaST2yLhuylNqW5otVPEySm9m7JZ0k6Qp3f7WkUyXtzd1vb7psobFC1+PdfV8arN5jZv9a94HpCtwmSZvWHe99FTA2U/pAhfZeaKJHW/mwlV9G8MIYQqtbJemO227TtW+6tKunuy03ZOkIZvZeSQ8queml7v4XFQ87StJ/lfTtkg5Iep+ZXSfpKyX3XZpJRgld7r4v/X2rmb1L0tmSbjGzk9NWrpMl3TpG2WIW4geqrVBavUIJXKGsDyxXFrTyCFwY05zrEnc/t8XD9kr6B3e/TZLM7GpJj1Eyzuu03P1Ok7Rv2ZMNPpDezO5jZvfL/lbSX7pL0lWSLkjvdoGkqtTZWigH0Dlrug1CuF5mSEIsE+6xKHDtvvCpRwQutifGwkk+tb1b0mYzW5sOqn+CpI+lw6FuN7Ot6VmLz1aN3DJGS9cDJb0rHax2lKS3ufvfmtmHJV1pZs+V9FlJP9Lli3a1c9X5djDWjjzVby5Dt+CFXhHN+ZtqiJa1bFUJfT/DPLAfJszsByX9jpJxW39tZte7+/e6+5fM7LWSPqyk+/Bqd//r9GE/I+lySccqGUC/cBC9NELocvdPSnp0yfIvSjqnyXPd8sU7BzsgN3mN/H2H2qGnfgAe6mzNmCogJtEdX9vABSAs7v4uSe+quO1PlHQnFpfvkLSpyeuYe7xj0U+yo/2HdbKk/lqgis+7yiDYOq9f9T7KHhvjQXbIUFN3/cQUtOqKcd+IUd3QRbciYvBGfea6RQPR+7Bm7To/6uHdjHO86/rLBi9/UyHN09W7i5+zuXZlV3bQyirYtsFrlQPhFA6iQx9oylqC5nKwowtyXAyWRxeKn9+51F9TNonQ1WXXX9VzFb/R1vmGS8UbljlWWIQvIC6LPqtjDF1Bt6IOXQ888Vhd/APhth4taxErC259BzUmaZwnwhcQtkW9K0W7L3zqoft3eZLYGy//TCfPhWpRh65VDBU+stcZKkyVvU7VB7e4nBA2fVOay21sDKJHV8rGDi+S/0Jf9nme8qXWYjfmBa9H06arsOvXXPV+VY9Z5b1xEJkHuiVWV+egCCxz8XM2Nw5cde6XPW/xB+ObXUtX3VafIV970f2KrU91nmOV97JscsdVzGkge+ho8WqPQIUu8Pmbp9mErlgrypDKPVRXKYbBOK/6uvwc8sVj3hZ93pruZ9TJ8ZlF6AopuEzBKnOV0doVHlq9gHaanE1Y98z4NvrsoUC3Jh+6CFwAVrX7wqc2qks40E1bWYAqO5tw2ZeZoYa1sD+GY9Khi8DVH1q7poXWrnBtQmIAAAwySURBVOWaBi/MU5fdh11ZNPUEhjXZ0EXl2L9VJogleIWH4AUs1/YzEuIxiVaw4U1yyogQd+65WnZaMwDEos30DhsvuyroY1JWvpDLOCWTC13sOOEheMWD1kfEro95qRbNp7WofuN4hKKouxePWb/hsP/ZwcO16NTmri9nAayKqzXEZ9Hgdqld/VL30jxM3YC6og5deQSuONQJXxIBbEyM7ToS9Uu46u6rq1wup8nE1gQvLBJ96KIyjNOyyokAhrFQp4zjD/zTkqSftvUL71cnZNU5W2/jZVfpYi0PSG0nLAXKRB+6EK+63wrpfsQQhjpYzv3M3WVddln46kPZtWnL6iCCE/pC6MKomjTHz/1ghf5wkO1PCBOELnv9YgsY0BdCF0a36ALfRbR69Y/xXOhC23FRYwi5bJgWQheCUlb5VQ26J3ghVlPef/NhizADHG5y83Rheqoqblpkusc6xSoIXMBihC5EYdXg1ceEiVMz53Uzxmn+U1rfxc8WgQsoZ+4+dhla27Jli+/YsWPsYmBgVQfIsu6aqgPbVLt22ppSAGir66DQJMjFuj/WnTwUcTCz69x9y5CvuWbtOj/q4d186bnr+ssGL39TjOlCdKrOeGwSHKY8pqYpAtf4YpyXjv0GaI7QhSh1MfPz3INXCAfNRa0iQ3b5hdQ6E8MZunQlAu0QuhAlLrXRXghhS1p+sJ77JVVCbP2KaRoIIESELkRnzgfiVcUSuPL3Y3uP3/rF2C2gG4QuzNpcuhhDCVtSd9eyaxvGyibjzX6HHiTGCF90JQLdIXQBExZS2OpaFgBWaQmLtSVtiPBFVyLQPebpAiZoTnOSrRoGYg4TfWzjsn0n5nUEhISWLmBC5hK0ipZdv7POoP1YdTHgfq77DTA0QheiEmNX0BDaHDRj7VpDtToBjIAFjIfQBUSs7QE0a9mp08LTdTDbfeFTe21Z6mKs16pCCLNdhauYWwGB0DCmC9EI4UAWgqGvI9nHQXeIbUlYWB3rEOgWLV0IXt8H6FimjRiz5aLYFblqa9JQB3NCQz2sJ2AYhC4EjdatbsfgrHJwLXtscVmd7cUBHsBcEboQrCEDV4itXU3DVghhps7g/L7HdAFAqBjTBaRCOaurzXitkELMGIPzASAGtHQBOWNe427VMxEBAGEjdAEluphwss1r1RVD0CorIy1cAObM3H3sMrS2ZcsW37Fjx9jFQA9CPjh3EcKWBa0YQhUWC3kfzmNfQ8bMrnP3LUO+5pq16/yoh3fzWbnr+ssGL39TtHQhSCHPll4MTHVDWJ0WLQ6A07Hs0kQA5ofQhWCFHLzymPkby8SyLwPoF2cvAsAACNUAaOlCUObaGsABeR5Ca/FivwOGRUsXMDIOfPMSwvbeeNlVQZQDmBtCFzAiDnzzNOZ2Z58DxkPoAkbCwW/extj+7HPAuAhdCMpcDgpzeZ9YbMhuPvY5YHwMpEdwQhts3BYHOdTV9z7PvgiEgZYuoAcc5NBUX/sM+yIQDkIXAEwUgQsIC6ELAALRZUgicAHhYUwXotDkADKF8WCYr7J9vek+TeACwkRLF4KUP2g0PYCMfcAZ+/UxPU32KfY/oDkz+3/N7F/NbKeZvcvMji/cfoaZ3WFmv5hbdp6ZfdzM9pjZi+q8Di1dmKRlB54+WsM42KFPVfsXLbtAJ94j6cXufreZ/aakF0v6ldztr5P0N9k/ZnYvSb8n6cmS9kr6sJld5e4fW/QihC7MUpen6BO2MCb2P2B17v53uX+3S3p69o+ZPU3SJyV9NXefsyXtcfdPpve5QtL5kqYbuq677rrbzOwzY5djROsk3TZ2IXp3uS26dfx1sLh8Qxh/HYyPdZBgPbAOpPbr4MFdF2QZv/OL777r+svWdfR0x5jZjtz/29x9W4vn+XFJfypJZnYfJS1eT5b0i7n7nCrpptz/eyU9btkTRx263P2kscswJjPb4e5bxi7HmFgHrAOJdZBhPbAOpLjWgbufN9Rrmdl7JT2o5KaXuvtfpPd5qaS7Jb01ve3lkl7n7neYHfYFu+zbti8rQ9ShCwAAoA53P3fR7WZ2gaTvl3SOu2cB6nGSnm5mr5Z0vKSDZvY1SddJOj338NMk7VtWBkIXAACYNTM7T0k34hPc/UC23N2/M3efSyTd4e6/a2ZHSTrTzB4i6XOSniHpR5e9DqErbm36qqeGdcA6kFgHGdYD60BiHbTxu5KOlvSetBtxu7v/dNWd07Mcny/p3ZLuJelSd79x2YvYPS1oAAAA6AuTowIAAAyA0AUAADAAQlckzOzTZnaDmV2fzUNiZvc3s/eY2SfS3yeMXc4umdmlZnarme3KLSt9z5b47fRyDDvN7DHjlbxbFevhEjP7XLo/XG9mT8nd9uJ0PXzczL53nFJ3y8xON7P3m9luM7vRzC5Ol89mf1iwDmazL5jZMWb2ITP7aLoOXp4uf4iZfTDdD/7UzL4pXX50+v+e9Pb1Y5a/KwvWw+Vm9qncvnBWunxyn4dYEbri8t3uflZu/pUXSXqfu58p6X3p/1NyuaTiHC5V7/n7JJ2Z/lwk6Q0DlXEIl+vI9SAlc8eclf5cLUlm9gglZ9E8Mn3M71tyuYrY3S3pF9x9o6Stkp6Xvtc57Q9V60Caz77wdUlPcvdHSzpL0nlmtlXSbypZB2dK+pKk56b3f66kL7n7BiWXcfnNEcrch6r1IEm/lNsXrk+XTfHzECVCV9zOl/Tm9O83S3raiGXpnLtfK+k/Cour3vP5kt7iie2Sjjezk4cpab8q1kOV8yVd4e5fd/dPSdqj5HIVUXP3m939I+nft0varWRG6NnsDwvWQZXJ7Qvp9rwj/ffe6Y9LepKkd6TLi/tBtn+8Q9I5VpjhMkYL1kOVyX0eYkXoiodL+jszu87MLkqXPdDdb5aSClnSA0Yr3XCq3nPZJRkWHZCm4PlpV8Glua7lya+HtIvo2yR9UDPdHwrrQJrRvmBm9zKz6yXdquQixf8u6cvufnd6l/z7PLQO0tv3Szpx2BL3o7ge3D3bF16R7guvM7Oj02WT3BdiROiKx+Pd/TFKmomfZ2bfNXaBAtPqkgwRe4OkhyrpWrhZ0m+lyye9HszsvpL+XNLPuftXFt21ZNkk1kPJOpjVvuDu33D3s5TMAH62pI1ld0t/T3IdSEeuBzPbJOnFkr5V0rdLur+SyT6lCa+H2BC6IuHu+9Lft0p6l5LK5pasiTj9fet4JRxM1XveqxaXZIiVu9+SVroHJf2h7uk2mux6MLN7Kwkbb3X3d6aLZ7U/lK2DOe4LkuTuX5Z0jZLxbcdbMkO4dPj7PLQO0tuPU/2u+ijk1sN5aRe0u/vXJV2mmewLMSF0RcDM7mNm98v+lvQ9knZJukrSBendLpD0F+OUcFBV7/kqSc9Oz9LZKml/1u00RYXxGD+oZH+QkvXwjPSsrYcoGTj7oaHL17V0HM6bJO1299fmbprN/lC1Dua0L5jZSWZ2fPr3sZLOVTK27f2Snp7erbgfZPvH0yX9fe6aetGqWA//mvsCYkrGteX3hUl9HmLFZYDi8EBJ70rHfx4l6W3u/rdm9mFJV5rZcyV9VtKPjFjGzpnZ2yU9UdI6M9sr6WWSXqXy93y1pKcoGSx8QNKFgxe4JxXr4Ynp6eAu6dOSfkqS3P1GM7tS0seUnO32PHf/xhjl7tjjJT1L0g3pOBZJeonmtT9UrYNnzmhfOFnSm9OzMNdIutLd/8rMPibpCjP7X5L+RUk4Vfr7j81sj5IWrmeMUegeVK2Hvzezk5R0J14vKbuMzRQ/D1HiMkAAAAADoHsRAABgAIQuAACAARC6AAAABkDoAgAAGAChCwAAYACELgCjmMDFlwGgEUIXgJWZ2ben13s7Jp3M98b0siTF+z3RzN5vZm+TdMMIRQWA0TBPF4BOpBNTHiPpWEl73f3/KbnPEyX9taRN7v6pYUsIAONiRnoAXfl1SR+W9DVJL1hwvw8RuADMEd2LALpyf0n3lXQ/JS1eVb46THEAICyELgBd2SbpVyW9VdJvjlwWAAgO3YsAVmZmz5Z0t7u/LT0r8Z/M7Enu/vdjlw0AQsFAegAAgAHQvQgAADAAuhcBdM7MHiXpjwuLv+7ujxujPAAQAroXAQAABkD3IgAAwAAIXQAAAAMgdAEAAAyA0AUAADAAQhcAAMAA/n+HjdV9urelaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 756x504 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "_ = psi.mean(\"t\").where(ds_xorca.fmask.isel(z_c=0)).compute().plot.contourf(size=7);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:xorca_env]",
   "language": "python",
   "name": "conda-env-xorca_env-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
}