GS_transport.ipynb 40.9 KB
Newer Older
Patrick Wagner's avatar
Patrick Wagner committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Once deleted, variables cannot be recovered. Proceed (y/[n])? y\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "%reset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load required modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load data\n",
    "dir=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/DATA/1_VIKING_K301/\"\n",
    "meshdir=\"/gfs1/work/shkpwagn/ARIANE/GRID/\"\n",
    "figdir=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/FIGURES/\"\n",
    "\n",
    "file_mesh=meshdir+\"VIKING20_mesh_mask.nc\"\n",
    "file_psi=dir+\"VIKING20-K301_1y_19490101_20091231_psi.nc\"\n",
    "\n",
    "data_psi=xr.open_dataset(file_psi)\n",
    "data_mesh=xr.open_dataset(file_mesh)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fix longitude\n",
    "\n",
    "def make_lon_monotonous(lon):\n",
    "    lon += (lon < 0) * 360\n",
    "    return lon\n",
    "nav_lon=make_lon_monotonous(xr.open_dataset(file_mesh)['nav_lon'])\n",
    "nav_lat=xr.open_dataset(file_mesh)['nav_lat']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Extract data and calculate area of grid boxes to weight spatial average\n",
    "\n",
    "xmin_nac=310\n",
    "xmax_nac=320\n",
    "ymin_nac=40\n",
    "ymax_nac=50\n",
    "\n",
    "xmin_gs=280\n",
    "xmax_gs=287\n",
    "ymin_gs=32\n",
    "ymax_gs=35\n",
    "\n",
    "def extractdata(xmin,xmax,ymin,ymax):\n",
    "    index=(nav_lat>=ymin) & (nav_lat<=ymax) & (nav_lon>=xmin) & (nav_lon<=xmax)\n",
    "    psi=data_psi['psi'].where(index,drop=True)\n",
    "    e1t=data_mesh['e1t'].where(index,drop=True).squeeze()\n",
    "    e2t=data_mesh['e2t'].where(index,drop=True).squeeze()\n",
    "    area=e1t*e2t\n",
    "    return psi,area\n",
    "\n",
    "psi_nac,area_nac=extractdata(xmin_nac,xmax_nac,ymin_nac,ymax_nac)\n",
    "psi_gs,area_gs=extractdata(xmin_gs,xmax_gs,ymin_gs,ymax_gs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average over latitudes and find maximum along longitude to get transport\n",
    "nac=((psi_nac*area_nac).sum('y')/area_nac.sum('y')).max('x')\n",
    "gs=((psi_gs*area_gs).sum('y')/area_gs.sum('y')).max('x')\n",
    "\n",
    "# PSI of 1992 is given in m3/s instead of Sv. My quick and dirty fix is to just convert this year into Sv.\n",
    "nac[43]=nac[43]/1e6\n",
    "gs[43]=gs[43]/1e6\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAADyCAYAAABAr8L8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsvXd4HOW5sH+/6r3LsoptFRv3XgDbGDC9JYBDhxASQgrpfIHwJb8TcnKSc8KXhOQkhMCBhCTUHHoxxQSDwTY2brjLRZZl9d679v398exIK2m2SbtarTz3de212tmZ2UfSzjzv05XWGgsLCwuL05eQQAtgYWFhYRFYLEVgYWFhcZpjKQILCwuL0xxLEVhYWFic5liKwMLCwuI0x1IEFhYWFqc5liKwsLCwOM2xFIGFhYXFaY6lCCwsLCxOcyxFYGFhYXGaExZoATwhLS1N5+bmBloMCwsLi6Bi586dtVrrdHf7BYUiyM3NZceOHYEWw8LCwiKoUEqd9GQ/yzVkYWFhcZpjKQILCwuL0xxLEVhYWFic5gRFjMCMnp4eSktL6ezsDLQofiEqKoqcnBzCw8MDLYqFhcUEJ2gVQWlpKfHx8eTm5qKUCrQ4PkVrTV1dHaWlpeTl5QVaHAsLiwlO0LqGOjs7SU1NnXBKAEApRWpqavBZOw0eJShYWFiMM4JWEQATUgkYBN3vVvg2/H4BVB0ItCQWFhZeEtSKINAopbjnnnv6X//617/mgQceGLTPwoULuemmm4Yd++tf/5pZs2Yxb948Fi5cyN///nd/i+tf9j4vz3XHAyuHhYWF11iKYBRERkby0ksvUVtba/r+oUOHsNlsbNq0iba2tv7tf/7zn9mwYQPbt29n//79bNq0Ca31WInte7rb4Mjb8nNrVWBlsbCw8BpLEYyCsLAw7rrrLh566CHT95955hluu+02Lr74Yl577bX+7b/85S/505/+REJCAgCJiYncfvvtYyKzXzj6LvS0y88tlYGVxcLCwmuCNmvIkZ+9foCD5c0+PeecrAR+etVct/vdfffdLFiwgHvvvXfYe88//zwbNmygsLCQP/7xj9x00020tLTQ0tJCQUGBT+UNKAdehthJ8nPraaQIqg/Dlj/AlQ9BWESgpbGwGDF+twiUUqFKqd1KqTfsr59USp1QSu2xPxb5WwZ/kpCQwBe/+EX++7//e9D2Tz/9lPT0dKZNm8YFF1zArl27aGhoQGsdfIFgV3S1wpF3Yc7nISETWk4j19DBV2DPU3BqW6AlsQgmulpg3wtgswVakn7GwiL4LnAISHDY9kOt9Qu++gBPVu7+5Hvf+x5Llizhjjvu6N/27LPPcvjwYYyuqc3Nzbz44ovceeedxMbGUlRURH5+foAk9iFH3obeDph7DTSWQEt5oCUaO2qPyHPRB5B3TkBFsQgidv4N3v0xtFTAym8HWhrAzxaBUioHuAJ43J+fE2hSUlK4/vrreeKJJwCw2Wz87//+L3v37qW4uJji4mJeffVVnn32WQDuv/9+7r77bpqbxZ3V3NzMY489FjD5XWLrc50JdOBliJsMU8+G+IzTyyJwVAQWFp5S/LE8v/czKN8dWFns+Ns19DvgXmCoDfQLpdRepdRDSqlIP8swJtxzzz392UObNm0iOzub7Ozs/vfXrFnDwYMHqaio4Bvf+Abnn38+y5cvZ968eZx77rnExMQESnTXfPCf8IclcHTD8Pe6WmT73KshJEQUQlsN9PWOvZxjjc0GtccgNALKd0FHY6AlsggGbH1wcgvM/hzEpsMLXxH3aoDxmyJQSl0JVGutdw55635gFrAcSAHuc3L8XUqpHUqpHTU1Nf4Sc1S0tg78AzMyMmhvb+eBBx7gvPPO45NPPhm0b2hoKBUVFWRmZqKU4t5776WwsJD9+/eze/dubr311rEW3z2dTbDtUfn5lW9Aa/Xg9wvfgr4ucQuBWARoUQbBwke/hX9cA96m7zaXDrjEtG1glWdh4Yqq/dDVJIrg2kehvgje/lGgpfKrRbAK+JxSqhh4DlirlHpKa12hhS7gr8AKs4O11o9prZdprZelp7sdsGPhDz59Arqa4ZpHZfX/yjcH3zAPvAzxWZBj/xfGTZbnYMkc0hp2/hWOvw+lXg4+MtxCC2+C8FjLPWThGcaCIXcV5K2B1d+H3f+QaymA+E0RaK3v11rnaK1zgRuB97XWtyqlMgGUpM5cDez3lwwWo6CnAz75ExSshYU3wsX/Acc2DFgInU1w7D1ZEYfYv0bxmfIcLHGC6oMS4AbY9Tfvjq09Ks8Z8+SithSBhScUfwwp+ZCQJa/P/7+QvRRe/y40ngqYWIEoKHtaKbUP2AekAf8RABks3LH7KXHxnGNvobH8TjjjUtjwb9JPqPAt6OsecAuB3TVE8FgEhevl+YxLYf9LYvV4Su0RiEqC2DTIPw/qjkJTqT+ktJgo2GwSH8hdPbAtNBzWPS6xg5e+GrD42pgoAq31B1rrK+0/r9Vaz9daz9Na36q1DnykxGIwfT2w+ffi8pm2SrYpBZ9/GKKTJMD12bOQOAVylg0cZxSVBYtFUPgWZC+DNT+EnjbY/6Lnx9YehbQz5O+Sf55sK/rQH1JaTBSq9kNnI+QOSTVOyYcrfgMlW+H9nwdENKvFhMVw9r0ATafEGnAsfotNg6sfgZpD4gqZ8/nB74dFQExqcFgELZVQthNmXiam+aQ5sMuLxn+1R0QRgBwbm265hyxcY8QHjMWVIwtvhKV3wObfwYFXxlYuLEVgMRSbDT5+CCbNhTMuGf7+9AvgrLvl5/lfGP5+3OSRWQSdzd5n7owGo0nezMtFmS35oiiGSg9CVp1N0lwvbbq8NqyCog/G9newCC5ObobkPEjMNn//sl9BznJJyqg+NKaiWYpgFFRVVXHzzTeTn5/P0qVLOfvss3n55Zdpb2/nlltuYf78+cybN4/Vq1cPSjUd1xS+CbWFcM4PBq/2Hbn45/C1TZC1ePh78RlSMekN3W3w29nw0a+9l3ekFL4FSdNg0mx5veAGqQnY/Q/3x9Yek2fDIgBRBG3VY34BWwQJNnuKsWN8YChhkXD93yEiFp67ZUxrUyxFMEK01lx99dWsWbOGoqIidu7cyXPPPUdpaSm///3vycjIYN++fezfv58nnngiOGYPay159cm5MOdq5/uFhELmQvP34iZ734q6uQK6W2HTb8Ym4NrdJqt3wxoAiEmB2VfBZ89Bj5vJcEbqqKMiyDtXni33kIUZ1Qfs8QEXigAkm+j6v0HjSXj562PWj8hSBCPk/fffJyIigq9//ev926ZNm8a3v/1tKioqBlUVz5w5k8jIICigPvYvqZJd9T0IHWEbqvgMUQTefIENxdHbAe89MLLP9YaiD6C3U+IDjiz5olysh153fXztEQgJE4VpkDQFUqdbiuB0ofBtyaDzlOLN8mwWHxjKtJVwyS/hyFuw6f+NTD4vmRBtqHnrR1C5z7fnnDwfLvsvp28fOHCAJUuWmL735S9/mYsvvpgXXniBCy64gNtvv50ZM2b4Vj5fU30IXrpTMhgW3Tzy88RNBlsvdNRLcNkTDEUw8wrY97+w4i6YYlpn6BsK10NkolxwjuSuEXfRrr/BguucH197RP5OoUOsvPzzYM+zknU19D2LicXe50QZXPTvnu1f/JEsHJKmeLb/irugbJe0eMleCjMuHLGonmBZBD7i7rvvZuHChSxfvpxFixZRVFTED3/4Q+rr61m+fDmHDo1j33HDSWmzEBoJt74kvsqRYtQSeDOgxmhJcel/iiJ56z7/mcS2PrmAZ1w0/GYdEgJLbpOL1lWjPSN1dCj550kaqqsq5d5uqD8hqaYHX5XXFsFHw0mxYLvb3e9rs0mgeJobt5AjSsFVv4MzvwbZ5gtOXzIxLAIXK3d/MXfuXF58cSDv/OGHH6a2tpZlyySvPi4ujmuvvZZrr72WkJAQ1q9fz+zZs8dcTre0VsM/rpYJY3e8BSl5ozvfoDYT8zyUoQpUqNQlXPgAvPJ12PdPSanzNWU7ob12uFvIYNEtsPGXUlB34U+Hv9/XI/1hZl0+/L3c1aBCxD007WzZ1tUqv8u+F6HhBDSXAw6ZRTc8DbOvHO1vFRzUF8Gb98B1T0JUYqClGR1GRXpHPUS4aRhZfRA6GtzHB4YSHi2ZRGOAZRGMkLVr19LZ2ckjjzzSv629XVYHmzdvpqGhAYDu7m4OHjzItGnTAiKnSzqb4KlrZfV+ywuQ4YO5Dv0WgRcB49YqiJskK/IFN0DWEokV+KMrY+F68e9Pd2JqJ2TBjIthzzNiPQyl4STYeswtguhkyaQq+kAsirfvh9/OgTe+LzeC/PPg3PukMO+m5+QY44ZyOnB4vfR1qvgs0JKMjq5WWUwAtNe53/+kPT6Q60F8IEBMDIsgACileOWVV/j+97/Pgw8+SHp6OrGxsfzqV7/i+PHjfOMb30Brjc1m44orrmDdunWBFnkwPR3wzI0ybvHm53znkx9J47nWalEEIMrg0v+Cv1wsxTVrf+IbuQwK35KAXXSS830W3ih1BsUfQ/65g98zyxhyJP88ybz6wxJROHOuHoh5OKbjag1h0dBcNprfJrio3CvPzV6mF483HJV3e737/Ys/gqSp8hinWIpgFGRmZvLcc8+ZvvfFL35xjKXxkvf/Q0rav/CE89XxSIiIkUCs1xZBxsDrqWfCvC/IPOAlX/TdBVR3HGoOSwWnK2ZcIh1F9784XBHU2ZvNpU43P3b+ddKMb+blsPRLED/ZfD+lxPpoPo0mulXYFYG3dSbjjUGKwI1FYLNJxpAzV+Q4wXINnY70dIoPfO41MM8Plkp8xsgtAoOLfgYo8df7iv5q4ktd7xcRIzGAQ69JTMCR2iPSU8mZRTFpthTbnfcj50rA4HRSBN3tUqgIE0ARnBz4uaPB9b41hyWO4G18YIyxFMHpSOF6yZdfcpt/zh/nxchKm82uCDIGb0/MgeVfgb3Pu87g8ZTWGvjkEciYPzj/3xnz1slFPrQuwFnG0EhIyD59FEH1QRngA8H/OzechLAo+dmda6hkqzwPTVUeZ/hdESilQpVSu5VSb9hf5ymltimljiqlnldKRfhbBr+itQSPgqnHzO6nJEMn71z3+46E+MmeWwQdDaD7hisCgFXflZTWDx8cnTy9XfD8rdBWC5//g2fHFKyVzBbHjqRaQ00hpPmoJiQhC1rKx6x61CcUfQjvjiBuYwSIk/MmhkWQnCsuUHeuoeZye0bc+I0PwNhYBN8FHJPofwU8pLWeATQAXxnpifV4uPl2t4nfuKvZp6cd8e9m65MbnzOaSiVzY9HN0irCHxgWgSe/g1FMFmsyhS5uklgF+/450N/HW7SWrJ1Tn8A1j5j3RzIjLBJmXQWH3hhoOdFeJ5aUzyyCLCm+C6bRnh//VmI33W3eHVe5V+Y3TDnTuxqT8UjjSYlbxaSI28cVbTVSWBkyvp0vfpVOKZUDXAE8bn+tgLXAC/Zd/oZMKfOaqKgo6urqAq8MeuwFJd5eGC7QWlNXV0dUVJT3B7/yDXjiouG+bYM9zwJ6dNXD7oifLMU2nihHQxGYWQQg7S7ComDTCK2CrX+EPU/DuT8aPETHE+ZdC90tMpkN3GcMeUuCvQ1JsGQOdTQOtFKuP+HdsRV7pVo/IVMsgmCygobSUCIV6DEp7i2CtlrzRc44w99ZQ78D7gXi7a9TgUattTGGpxRw0pPVNTk5OZSWlhLwwfbt9dIwLawF4pp8dtqoqChycnIGNvT1SrGSu5VF+W65YX3yCKz6zuD3bDbY85TMSvXETz5SjBTSlkr3hUOt1fZjnCiCuHSZjrb1jzJAxhu3zJF34N3/T1I4z73P8+MM8s6V+Qr7X5SGdP2KwIeuIRD3wRhUj46aoxvEggEpjpvsYcFgX6/ECJbfKTOubb2Shz80QSAY6GiQ4fPJ08QyML6/zjAsgnGO3xSBUupKoFprvVMpdZ6x2WRX0yW9Uuou4C6AqVOH+9fCw8PJyxtlFawvePRcqNgD0Slwb5Hz1s2j5a+Xill9yS+c72Prg4ZiURgf/KesgB17m5zcLO+f/2P/yGjg2GYifabrffstAhc3hZXfgU8fl1jBuv/xTIbqQzJJLXOBDNMZiWkeaq8D2POMxIFqj4p1kuhhvxh39FsEQRI8PfyGFM11NEiVsKfUHpEmf5MXDFThNpcHpyJosGcMJU2Ta776sOv922ogeZnrfcYB/nQNrQI+p5QqBp5DXEK/A5KUUoYCygFMrwKt9WNa62Va62Xp6ePUtOrrlfSw6GTxFfqrSrS7XfrXlH7qer/mMpkjvPoH8vrtHw1+f/dTEuCafZV/5DToLyrzIHOotUoKqyLjne8Tlw4rvgr7X4CaI57JsP6HUqJ/47PuWwC4Yt46cXMdeVtuaKkzfOfvjUmFkPDgcA31dkl9xJzPi9zeKAIjUJy5QCwCCN44gXGNJ02Vv4PbGEFwuIb8pgi01vdrrXO01rnAjcD7WutbgI2AMdrqduBVf8ngd+qPy0pn/vXyuny3fz6n5jCgoc5NwNS4OPPWiCvk8BtSSQvSTuLgqzB/ndwg/Yk3jeeMGgJ3ltTK74jC8CRWYOuTnkLzrnU+DcpTpp4N8Zky3L72iO/cQiAKJSEzOCyCE5vEBTrrSsn88SZGULlX/nepM+T3BcmWCkaMGoLkaRCTLH8TZ8kZPR0SYwoC11AgQtn3AT9QSh1DYgZPBEAG31BlH2u44AZZ2flLEVQflOf2Otd5y8bFmZIPZ98N6bNh/b0SyN7/oqxsF9/qHxkdiUyQC98Ti6DNpJjMjNg0sQr2vSApnK6oOy5B/MkLPJPXFSEh4mI7tkFWg74KFBskZAdHOuXhNyAiThYZKfneKYKKvdLHKjRMivFUiH/bTJzcMrJxqZ7QcFKs6uhkcQ2B82uyzd6PKMZSBABorT/QWl9p/7lIa71Caz1da32d1tpFruM4p+qA9JOZPE++6P5SBFUHB352ZZLXF0nefUK2tFi+8rfQVCLDLXY/JUPWs8YgKKmUZA55bBE4CRQPZeV3xJrZ/pjr/YyeNpk+UAQg7qG+bimI8qVFAPbq4nHuGrLZxLKcfqGk1abkQXOp6zRlA61lVojxvzCUgb8sgu42+PvnJWXYHzSehGR7zDImVZ6duYeMtODT2TV0WlB1QFaIYZGSn16+xz+FZdUHJQcbXLuH6ovEZDV82NNWwqJbYfN/i6tk8a3+C2YPJd7DkZVG51FPiE2FqWdByTbX+1V8JvOH02d5dl53ZC8d6HfkF0VQPr4LEst2yv9p1hXyOiVflKInMbGGYsmycbTOEjL9ZxGc2i5Ku3C99ymuntBoTx0FSR8F5ymkhkVgKYIJTtWBgdbNWYvlC+9NEM1Tqg/KIBUV6kYRnJCL1JGL/h2iEsRyWXCD72VzRlyGe4ugr0cuIk8tAoCc5TL/1VWL6sq9Yv34akqYUvK3C4t23mxupCRkS5zJXc+aQHL4Dfn+zLhIXifbs/U8udGaWWfxWf4LFp/cLNdJSKh7y9FbtB6iCOwWgVPXkGERWK6hiUtHAzSdGqwIwPfuobY6WY1lLpTVvjNFoLXkdg9VBLGp8IW/wJUPje0X0hOLwLhQvEkjzF4mq9GKPebvay0+aV+5hQzW3Avf3AoRsb49b38twTh2Dx1+U5qmRSfLa+M75smip+IzuTFPcph1ET/Zf66h4s1yrcy9Fnb9Azp9WPHfViOxp2S7IjBiBJZr6DTG8NtnzJfnSbPFP+9rRVB9wH7+ObIadaYIWqvkSzpUEYD0zVkyxm2x4zKkstjVKD93VcVmZC+VZ2eptE2lcmH6IlDsSFjE6Ke3mTHeawlqjkgLlZlXDGyLTZPAcYMHFkHFXnHRhTtUySdkykKqp8O3svZ0QNkOGQBz1tclY2fPM747f3/qqKeuoRqxIn29ePADliIYKVX2G7RhEYSGSwl9uZOV6og/x65w+hXBcXN/srE6Sx4HRXYw0ILZVfM5d1XFZsSmirJzNhe43xWx0PNzBpLxbhEUvinPjqM5lRKl6IlFUGlinfXXEvg4TlC6Q+ID01bLgmHKmbDtz+aT5kZCQ7E8G/GisEhRiO1O3HpGDcFYxeVGgaUIRkrVfjENHXvOZy0Wl4Uv+6hUHxSTPH4ypBbIqt/sAupPHR0nisC4ubtK43PVcM4V2cvkojdTiBV7AeWbsZtjQX865RhaBDWFUPi2Z26Tw2+KUk3MGbw9Jd+9Imipkv/xUOusv5bAx3GCk5sBJQkFAGd9Q6yWI+/45vxGDYHjoKRoF/2GgqS9BFiKYORU7ZebjaO2z1osBSbuCr+8ofqg+FeVGghUmp2/vkh8seNlHJ5XFoGXrQZylst5zVbRlXslkysIzHFA0injJrtWBFUHfOtGefFOePYG+FUuPHExbPxPKPkEersH79dSKQp31pXDz5GcZ5/f7GK17SyNN96hx5IvKf5YrHJjaNCsqyAhBz75k/n+tj7pDuApDSelJiAybmCbqw6k7cFRVQyWIhgZtj7pZTN5/uDtvg4Y22zyORlz5LU7RZA01XeZMqOlv/GcK4ugWopzvK10znERJ/BHoNjfuKolaK+XflZb/uibz2qrkxv0oltg9fekAdymB+Evl8AvMuA3s+Hxi+B/vwSvfxfQMnZzKCn5YOuRmIwzjNYSQ68TY5HgS9dQb5d8HxwngYWGSRFi8UdQuX/w/ofXw2/nwAtuxpY60lgyECg2cNWBNEjaS4ClCEZGQ7G4aIa6H9LOgPAY3ymCphKxMCbZFUF8lgSfzCZ21ReNH7cQyAUSEu76YvemhsCRjPkSmB8aJ2ivl0InXweK/Y2rkZXlu+WGW/yRbz7LOM/SL8EF/wZffV+aJV7/d+nuWnC+KOaKvXB8o/wtzdxsxnfNVcC4cq9YDkM70EYlynXiy1qCsl2Shjtt1eDtS74o18y2R+R1W51YRM/dJNfwoddkseUJxhwCR2JSzdNHtQ4q15A1vH4kGK0lhl4goWFy4XiqCI69B+//QiZxzTUZy9CfmWT/nJAQiRMMtQi0lhhBzjjqcqiUxAlcpZB6U1XsSFiE+K2HKgLH5mbBREK2DAsyo3yXPJd+KnUXo7X4TnwIEfGDK8yjk6WZ3JzPD95Xa3mYBTsdU0jzzzP/rIrPzIP2Skn/Jl+mkJ60z0kYOhIyJgUW3QS7n4acFfD+zyVj6bz7Yekd8N+LZNDO1U7cRwa2Pmg8BbM/N3h7tBPXUFezBK4ti2ACU3VAAnxmlatZi2Ul5Mr32FwhpvdT6+RC3/qw+X5GjyHHzzFTBEaPdLPU0UAS76aobKQWAUicoGLP4AE8/a6IYFMEWWL5mQVvy+yKoKd9wOc+Gk5skptlqAdrQKWcd1qNzxKrzFlRWWeTWM7O/hcJPi4qK94ssTQjpdORM78OfV3w+ndE6d71IZz3I/l+Lr4V9v7TvXXSUiGWmZlrqLNp+PUeRFXFYCmCkVG5X/z1Zr7trMVy0daatEu29cG2x+DhFeKjPP8nMjmrdLt5uX71QZl1GpUwsC11ulxgjl+88ZY6ahCf6R+LACRO0Ns5YJ2B3CgTp5jfDMYzCS6Cp2W7BmZLl3wyus9pKpNFRN6a0Z0HREEk5zrPHKrcJ8/O0njjXQTI22rhD8vENeUJfT3SWiJ3lfn76TNh7U/gop/Dnf8aPFDnrG/KzOxtf3b9GUNrCAz6+w0NSSENoqpisBTByKjaDxlOpjM5Cxg3lckIybd+KDnO39wK5/4QFtrbPhx42eRzDg4Eig1Sp0uAz0hlg4GLcbxZBK7aTHS3S8FP3AhXTNl2N5ije6hib/DUDzjibGRlc7lkR828XG66JVtH9zknNslz/rmjO49BSt5Abv1QTtplzVxk/n58pnw3zFKAT3woRWyvf8d1QaJB+R7oaRseH3BkzQ9lYt9QSyglT9w9O/4KXS3Oj3ccSOOIUW09NGAcRFXF4EdFoJSKUkptV0p9ppQ6oJT6mX37k0qpE0qpPfaHk2/KOKWzWW7CzvLUU6dLkYmjIqj4DB6/QKo01z0Bt70sLh6Qm3fWEmkT7Uhvt1wMk0wUAQwOGNefAJR/x0+OhPjJ4j8dmpYI0n4aRm4RJE2VHHxDEXTZ03aDzS0EA3n1Q1fIhlsoe6nMRSj5ZHTN6U5sEp/2JB/VWBjtqM1kKnxT3HfOFH1ClrhrzAKtpz6V3kaNJfDhr9zL0R8fcKEIXLHqO+Ja3fV35/s0ngTU4Il/4LwDqaUI+ukC1mqtFwKLgEuVUvZKD36otV5kf/i4FNfPGBkGziyCkBBZBRmK4Mi78JfLJMf/K+/A/C8MD77NWyfKotbB9197RFb+QxWOWQppfZGsKh3L+McDxk3ezD00kqpiR5SS4HiZXRFUHQB08AWKQVbHMFwRlO+ytzmfL0VSbTUjb2qotSiCvHN8N2EtOU9W4kPn9jaVyvffLO3UwPidzbLKTm2DKWfB4ttkVrVRxe+M4s2QNnMU1uVSqUbe+qfBMSdHGk6KzGGRg7c7azPRP4sgdWQyjTH+nFCmtdZGi8hw+2Mc99r1EGcZQ45kLRIf6Sd/lsKdtOlw53vOj5l7jTwfeGlgm6FwJs0evG9MipijQxXBeEodNYh3MbLSk1nF7shZJn+H9vqBQGowWgRhkbJyHOoaKtslFmF4lFgEIENXRkJ9kaTW5vnILQTOm88ZU/HMCtEMEpy0mejpkP/llOX2zrmJUs/grFq/r1csJWfxAU9Z+W35+5i5aMG8hgCcdyBtq5HW8WERo5NrjPBrjEApFaqU2gNUAxu01kYj+V8opfYqpR5SSkU6OfYupdQOpdSOmpoaf4rpHVX75cs5tOTekazFYva+fR/MuAS+tH7A/DcjMRumrpRxiAbV9qE3qSb971OGZA41nBifisDwfZvVPYyk4dxQjDhB2S7JIIpJHbjBBBtDawm0Fosg257mmXaGuHVcBYzfvEfSJM048aE8+1QRGO2ohyiCw2+K5ZruYpqbsUgYZgXtFkt4ypmy6Lnkl5I6u+tJ8/NU7pVY00jdQgYzLharYst/m7u6zGoIwHkH0raaoHELgZ8Vgda6T2u9CBlSv0IpNQ+4H5gFLAdSkNGVZseOz+FxxZAlAAAgAElEQVT1VQfELeSqkdTUsyVXe8XX4ManB5ekO2PetVBzaKB2oOqgfeiNyYrCaD4HErNoqxl/gWIQayYqaeAm5EhrNaBGN8Yve4mco/TTgUBxEDT4MiUhe/BNsb5I0hKNfH9l76HjLGBcvhs+fRzeug9aTRZOJzZJyqcRm/IFSVPF5elYVNbRKEVrxhAbZ8Q5qS4+ZV8r5qyQ5wU3SJbThgfMq9RPbpZnx4rikRASAiu/JZZ80QeD3+vrEWttaKAYICIGwqLMXUOWIhiM1roR+AC4VGtdYXcbdQF/BVaMhQw+wWazZ/K4CbYlZsN9xXD5gzIgwxPmXC21CUbQuPrg8ECxQep0MWO72wcuwvGoCEJCpdjo+PvDV1mtVZJa50k+uzMi40XZlGyxt/wIQreQwdA2E46BYoOpZ0P98eE+eYDtj0sFbW8HfPDLwe/ZbHDiI7mh+lJRhoZL8NTRIji6QVb0rtxCIAuc2HQTRbBdLN5Yu8tFKbjiIfm93vm/w89TvFn2d2z+OFIW3CAJCO/8WIrHDJpKZQaGmWsI7NXFJumjQZI6Ch4oAiXMV0pdopRao5TyKPqhlEpXSiXZf44GLgQOK6UyjfMCVwP7nZ9lnNFUImaoJ50tvb3BxaXLhbr/RVkJNp0anjpqYKzq6ovGbw2BQcFaudhrDg/e3lotF91oyV4qNzlbT3AGig0SsqCzUWbugriFwqIHFxMacYKh7qH2etj/Aiy8EZbfCTufHNw2ofqgNEDzVdqoI8l5g4vKCt+U/2u2B1Xu8UNGVmotimDKmYP3S5sO5/wf+R0fPgv+tHLgcfxfo48PGIRFwucflnjAo2vg2L9ke6OT1FEDsw6kE8U1pJTKVUr9CTgO/A64A/gBsEkptVkpdZv9Zu6MTGCjUmov8CkSI3gDeFoptQ/YB6QB/+Gj38X/GLMGMua73m+kzFsnK/w9z8prVxYBSJxgvLWfHkrB+fI8tIXCaKqKHclZTn8OwuQgrCEw6K8lsN8Yy3aJq8txQZG5UNwQQxXB7n9Icd2Kr8K594ml9O5PBt436gdyz/G93I7tqHu7xCKYeZlnmUlD20zUF4nCmmLiJFj9PSn+Si2Q77rxmHmZKD9fccbFcNcHYmE8tQ4+fHDgGnPW2XdoB9K+XlHOQWQRuFq2Pgg8AnxLaz0oZG9f1d8C3A48aXaw1novsNhk+9qRChtwTmySGgF/rTxnXQlv/AA++o28dqYIDDdQ3TEp6ImdJBf/eCRpqgS8j2+Es+8e2N5aYx4I9xajv1JE3Ph0j3lKfxZNudSDVHwGy4Z0xgyLkJV2iUPmkK1PYgPTVg1YqmvuhXd/DEffgxkXyvc2JX94DrwvSMkTS6a9Xobcd7e6jw8YJGTKMQantsuzmSIIi4RL/3P08nqCkeX3xvdh4y8gMkFiIYayHkpMykAlNdiVgp4YFoHW+nqt9cahSsD+XoXW+tda6yf9Kt14o79Pi59aPcekwPQLpNgqIt75CiQyTgJ/dcfNB9aPNwrWSq/43i55rbXvLIL0WaIEJs/3XX58IHDs0V9zSHzijo3hDKaeJYHxLntm9tEN4spY8dWBfVZ8VVw27/4YejoloOrLbCFHjO9ewwnJFgqP9fyz4rPEAjC+F6Xb5aZr1sNrrImIhWsehSt+IymtydOcu3ujUwanjwZZewnwLEawSyl1r1Iq1//ijGOay6XS118XlMG8dfI8abbrwJ7RfG681hA4UrBWbmyGS6OzSdJrR5M6ahASChc+ACu/M/pzBZL+6uIyh0CxmSI4W3rjGIV02x8TF4tjcDYsUnLwaw7Dmz+QTpi+6C9khhGbqiuCwvVigXha2Gj8zkYq8antYuF5mmDhb5QSt9PXP4Lr/uZ8v5hUsYqMIT1BVlUMnmUNXYe4kF5VSm1VSn1PKeXERprAGH5Wf11QBjMvk17t7txPqdPlQm8pH/8WQe5qmU1gxAlGW1U8lBVfHTxTNxiJiJVU2+ZyCRRHJZn/X6csB5Qo1brjEixdesdwK3X2VVKbssdeV+A3RZArz/tfkBv6TA/dQuBQUV0hadBVBwbSRscTk2a7vh5jUiSrqLNJXgdZ51HwQBForY9rrX9pbxXxZWApcNLNYROPog/FBHTWWsJXRMaLf/L8H7veL3W6rPRg/CuCyDjJBOlXBD6oKp6IGLUEZbukKNHMIoxKlO9gyVaJDYSEy5CZoSgFl/xCfs6Y5z83RUSM3NCPvC1+9DMu9vzY/jYT5fZYgTaPD4x3hlYXT1CLAKVUjlLqB0hgeA7g5i41wfBHnxZXZDjpq+6IkTkE4981BJI9VLlXgsSWIjAnwR73qTpg7hYymHa2NGbb/TTM+Zz01Tcje4lU5p57r3/kNTAWIrmrBrpxekJ/++0Ke6BYja/hSp5iVBcbKaRtNaIUo5ICJ5OXeBIj2AysB2KB27TWS7XWHrQEnED092nxs1vIGxwVwXitIXCkwJ4sVvSB711DE4WELKgtlBiAYyHZUKaeJc3euppgxV2uz3n23cMnj/kaYyHirohsKNHJMtympUICxZNmDx9rGQzEDGkzYRSTBVHygidVT1/TWgdP0Zc/6O/Tcl5AxRhE8jRZdUTGB8cglsyFsnI6/r6sYEPCg2rFNCY4pieaZQwZTLE38Z08f3jxVSBImylV8TMv8+44pSRg3FwmFs68a/wjn7/p70BqKILgai8BLhSBUupyYL+hBJRS/xdYh8QHvq+1Pn3iBP7o0zJaQsMlUBcsKyjHdhMFa8UtFEQrpjHBcJXEZ7pvUnjW3TDjovHRW2n5V8Rt6izd2RXxmVIZ3tU0PpTaSDBzDQVR6ii4tgj+E1gJoJS6AgkU34IUiT0KXOp36cYDNpsogunj5KJz5Nz7xt8MAlcUrJVW2yc2WfEBM4ybvytrwODSX7rfZ6yIiB2YzOct8ZkDjfTGY8aQJ0TGi4Xr6Boab0Oi3OBKEWittb3xCdcCj9vbSG9TSn3N/6KNE6oPiqb3R5+W0WKMuQwWjHYTzaWD58ZaCIZrKHuEN9VgxLCColPGl8XtDUqJe6jfIgg+15Ar2zxEKRVj7yd0AeDYLMZ0hsCExJ99Wk43EnPEnwxBd6GMCWkz4dwfwaJbAy3J2GGkkE45c/xZ3N4Qkyoxgu52abMxgVxDfwB2A03AUa31dgCl1ELAyUTyCciJD/3Xp+V0pGCtZMZYGUPDCQmB8+8PtBRji9E+esrywMoxWow2E+3BV0wGrnsN/Q9wEXA3g+MBtUi8YOLT1yv9zv3dVuJ0wkgjtRSBBUjNTEjYwPciWDE6kAZhMRm4zhqaorUuAUoct2uty+zvKyBTa13u5PgoYBPiRgoDXtBa/1QplQc8h0wn24XUJnT74pfxOeW7Zf7AeKofCHby1kgl7IyLAi2JxXhg0my4vxTCowMtyeiIsVsEQdheAlzHCH6vlHpeKXWzUmqmUipFKZVlH07zU+BjwFVj/i5grb01xSLgUqXUWcCvgIe01jOABuArPvpdfE9//YClCHxGeBRc9fvgqIa2GBuCXQmAuIY66geKJYMsRuDKNXQt8AtgIfAEMlzmHeBbSC3BhVrrd1wcr7XW9l65hNsfGlgLvGDf/jdkStnY0VI1MHnIHSc+9G+fFgsLi4lBTKqM6DRGxwaZReCystg+XGbvSE+ulAoFdgLTgYeRaWeNWute+y6lgGknU6XUXcBdAFOnjqBQxRmv3g3HNsC3dsoACmf0dELJNimWsbCwsHCFUV1cUyjdgyNiAyuPl/i1tFNr3ae1XgTkIEPqZ5vt5uTYx7TWy7TWy9LTfaRdiz8WJQCw/VHX+5ZslZ75VqDYwsLCHUYH0prCoPQgjEmNv9a6EfgAOAtIUkoZlkgOYBps9oMQ8N7PpFXE3Gukc2NHo/P9P/mT9MLJXT0m4llYWAQxRpuJ+qKgcwuBHxWBUipdKZVk/zkauBA4BGwEvmDf7XbgVX/JMIjCt6TD4Xn3wervS/fGXX8337dkGxx9VwZmR8aNiXgWFhZBjOEa0n0QMwEtAqXUu55sMyET2KiU2osEmjdord8A7gN+oJQ6BqQigWj/YuuDf/27tG5edKt0wpy2Wsb89fUO33/jf4hWd9fi18LCwgIGdwAOQovAVR1BBBAFZCil4gGj/jsBcBu9tQeahzVN0VoXIfGCsWPvP2Ug+HVPDgygPvub8NzNcPh1cRUZFH0obSUu/a+gC/hYWFgEiMhEacWtbRMuRnA3cACYZX82Hu8Af/a/aD6itws2/lKsgNkOAzrOuFQ6BH7yyMA2rWHjLySOsPSOMRfVwsIiSAkJGYgTBKFF4KqO4CEgF/ip1nqq1nqK/TFXa/27MZNwtOx8EppK4IKfDu5/HxIKZ34dTm2D0p2y7dh78vrcHwZXe2cLC4vAEzMBFQFI+idjXfDlS7pa4MMHpXOoWS+TxbdCZIJkCGkN7/8ckqadXt0fLSwsfIORQjrBXEMGG5RSfh566ic+eUS6AV74gHmL28h4WHwbHHwFtv8PVHwG5/0IwiLGWlILC4tgZyK6hhz4FvCyUqpDKVWvlGpQStX7WzCfkDodVnwNcpY53+fMuyTA89a9kDoD5l8/dvJZWJxG1LV28cBrB+jo7gu0KP4hJlmeJ6giSEP6BMUB6fbXwfGbzrsWLn/Q9T7JuTDrCkBLL/hQl103LCwsRsj7h6t5cksx7x2qCrQo/iE2XTKHDBdREOH2rqe17rMPsjdacH6gtX7bv2KNMRf+DCYvgDnXuN/XwsJiRJQ1dgCwsbCaqxZmBVgaP7D8Tpk3HYSuZbeKQCn1C2AV8Ix9071KqdVa65/4VbKxJLUAzr030FJYWExoyhpEEXxYWIPNpgkJCeLRlGYk5sgjCPHED3IVsNieQYRS6i/IQJmJowgsLCz8jmER1LV1s6+siYVTkgIskYWBp72GEhx+jveHIBYWFhObssYOVk1PRSlxD1mMHzxRBA8Cu5RSjyulngB2IFPGLCwsLDzCZtNUNHYyLzuRRVOS2FhYE2iRLBxwqwi01k8Bq4H1wJvAGq310/4WzMLCYuJQ09pFd5+NnKRozp85ib2ljdS1dgVaLAs7nrqGliKzBM4ElnhygFJqilJqo1LqkFLqgFLqu/btDyilypRSe+yPy0cmuoWFRbBQag8UZyeLItAaNh21rILxgidtqP8AfBc4ChwDvmPf5o5e4B6t9WxEidytlJpjf+8hrfUi+2P9CGW3sLAIEoxAcXZSDHOzEkiLi2TjYUsRjBc8yRpaC8zTWmvozxpyO8dYa10BVNh/blFKHcLJfGILC4uJTZmDRRASojhvZjobDlbRZ9OETrQ00iDEE9fQEWSkpEEmsN+bD1FK5SKzCbbZN31LKbVXKfUXpVSyN+eysLAIPsoa20mMDicuUtae58+cRFNHD3tONQRYMgvwTBEkAoeUUu8ppd5Dxk0mKaVeUkq95O5gpVQc8CLwPa11M/AIUAAsQiyG3zg57i6l1A6l1I6aGsuEtLAIZsoaOshOiu5/vXpGGqEhynIPjRM8cQ39YqQnV0qFI0rgaa31SwBa6yqH9/8HeMPsWK31Y8BjAMuWLdMjlcHCwiLwlDV2MC11YOJfYnQ4S6cms7Gwmv9zycwASmYBnlkEW4D3tdb/AoqBSOBDrfW/7NtMUUopZB7xIa31bx22Zzrsdg1eupksLCyCC631MIsA4LxZ6Rwob6aquTNAklkYeKIIPgKi7TfwD4FvAH/x4LhVwG3A2iGpog8qpfbZh9qfD3x/hLJbWFgEAU0dPbR195GTPFgRnD9zEiC9hywCiyeuoRCtdbtS6svAH7XW/6WU2uPuIK31xwwMvHfEShe1sDiN6K8hGGIRzJocz+SEKDYWVnP98imBEM3CjicWQYhSajlwMwP+/FD/iWRhYTGR6K8hGGIRKKU4f1Y6Hx2tpafPFgjRLOx4ogh+APwMeFNrvV8plY+4iywsgp7Onj6a2nsCLcaEpsyJRQBw7hmTaO3q5bNTjWMtloUDnvQael9rfbnW+hf210Va62/6XzQLC//z4NuFXPGHj7DZrMS0kfLhkRq+8MgWunrNR1CWNXYQFR5CSuzwgS3zcxIBKKxq8auMFq7xpMXEdKXUn5RS65VS7xqPsRDOwsLfFNW2UtrQwW6rsGnEvHewih0nGzhY3mz6vpExJImEg8lMiCI6PJTj1W3+FtPCBZ4Ei19A0kCfAibo1GmL05VaewfMt/dXsnRaSoClCU6O2Ffzu0saWTx1eKOAssYOspNjTI8NCVEUTIrleE2rX2V8ZXcZU1JiWDrNamRghicxApvW+g9a6y1a623Gw++SWViMATUtogjeOVCFvZ2WhZccrZab+G4nfv6yxuE1BI4UpMdxrNq/iuBnrx/g8Y+K/PoZwYwniuBVe7uHdKVUgvHwu2QWFn7GZtPUtXYzKT6Skvp2DlVYfmpvqW3tor6tG6Vgd8lw91p7dy/1bd3DaggcKUiPo6yxg45u/zgc2rt7aWjvodIqXHOKJ4rgTuD/Q+YUH7A/rGpgi6CnsaOHXpvm+mVTUArePlAZaJGCDsMttGZGOqUNHVS3DL7Zljc6zxgymD4pDsBv7iFDhqomSxE4w5OsoSkmj6ljIZyFhT8x4gMzJ8ezfFoK71qKwGuOVsnN+wZ7QdjuksHuIceBNM4oSPevIihrFAVQ3dJlZYc5waMJZUqpWUqpa5VSNxsPfwtmYeFvjPhAWlwkl8ybzOHKFk7UWtkr3nCkqoWEqDDWzppEeKgapgjKPLAIctNiCFFwvMY/f3ujjqHXpqlts8ZjmuFJ+uhPkC6gfwYuA34HfMHPcllY+B3DIkiPj+SSuRkAvGNZBV5xtKqVMzLiiQoPZU5W4rA4QVlDB2EhioyEKKfniAwLZWpKDMf9FDA2XEMAVU2WIjDDE4vgBqQ5XIXW+jZgIZ6lnQYNGw5WccFvPqCpw6owPZ0wLIL0uEhykmOYn51oKQIv0FpTWNXCjIx4ABZPSWJvaRO9Du0iyho7mJwY5XYKWUF6nB9dQwOKwAoYm+OJIujQWvcBvUqpeKASyPevWGPLx0drOF7TxvOflgRaFIsxpKa1i4jQEBKiZV1zydwMdpc0UmkFFT2ipqWLpo4ezsgQH//iqUl09PRxuHIg+8qs/bQZBZPiKKpto88PPvyyxg5yU6WOwVIE5niiCHYrpZKQ1tM7gO1IBpFLlFJTlFIblVKHlFIHlFLftW9PUUptUEodtT8HvMLjmH0l8uTm4kGrGYuJTU1LF2lxEf0Vr5fOmwzAuwctq8ATjtgDxTPtFsESezGZYz2BFJO5VwTT0+Po7rVR2tDucznLGjpYkJNEaIiyMoec4FIR2IfLPKC1btRaPwxcAXxNa/1FD87dC9yjtZ4NnAXcrZSaA/wI+JfWegbwL/vrgHKsupXMxCjKmzp5a791E5gI/HPHKV7ZXeZyn9rWbtLjI/tfT58UT0F6rOUe8hAjddRwDeUkR5MWF9EfJ+jps1HV3EmORxaBTC/ztXuoz6apbO5kSko06XGRlkXgBJeKQEup5RsOr49prd1aA/Z9K4x9tdYtyKzjbODzwN/su/0NuHoEcvuM5s4eqpq7uPWsaeSmxvDExycCKY6Fj/jj+8f4y2bX/0uxCCIHbbtk7mQ+Kaqnoa3bn+JNCI5Wt5AcE05anDSTU0qxeGoye+yZQ5VNndi069RRAyOF1NcVxlXNnfTZNFlJ0WQkRo1oGtp7B6sm/H3BE9fQdqXUktF8iFIqF1gMbAMytNYVIMoCmDSac4+WInvK2oxJcXx5dR57TjWy86TVgCyYae/upaS+fVC2iBm1rV2DLAIQ91CfTfPeoSonR1kYHKlqZUZG/KBmcounJlFU20ZDW7fDQBrzPkOOJMVEkBYX4fPmc44FbZMTIkcU/3nsoyJ+/sbBCX1fcKoIlFJGZtBqRBkUKqV2KaV2K6U8sgrs54lDBth/T2tt3p7Q/Li7lFI7lFI7amr8N8rOWIEUTIrjC0tzSIwO54mPrZ4kwYzxP61t7aazx7xtQZ9NU9c63CKYn51IVmIU7xywFIErtNYcqWrpDxQbLJ4icYI9pxqdDqRxRr4fMofKBimCKK9dQ8bvCfDAawcmbEGaK4tgu/35amAmcDlwHVJDcJ0nJ1dKhSNK4Gmt9Uv2zVXGAHv7c7XZsVrrx7TWy7TWy9LT0z35uBFxrLqV8FDFtJQYYiLCuGnFVN7eX8mpet8HrXzN/rKm/i+pxQCFDlkrFU5WgA3t3dg0wywCpRQXz53MR0drnCoRC6hq7qKls5cz7PEBgwU5iYTY+w4ZhVyZic5rCBwpSI/jWE2rT5v/GYrAcA21dPbS3t3r8fE1LV00tvewPDeZfWVN/HPHKZ/JNp5wpQgUgNb6uNnD3YntgeYngENa6986vPUacLv959uBV0cou084Vt1KbmosYaHyp7h95TRClOLJLcWBFMstbV29fPEv2/n31w8GWpRxh6NyNG5GQ3GsKh7KOTPS6Oq1scukiZqF0B8onjRYEcRGhjFrcgK7TzVS1thOenwkUeGeTbadPimOxvYe6n0Ynylr6CApJpzYyDAm24vavHEPGamw37/wDFbkpvDgO4UTst7IlSJIV0r9wNnDg3OvAm4D1iql9tgflwP/BVyklDoKXGR/HTCO17T2N70CyEyM5ooFmTz/6SlaOsfvP/zpbSepb+u2siBMOFLVSlJMOIDTOIFjVfFQVuSlEBqi2Hq8zn9CBjmGIhjqGgKJE+wpaeRUvWc1BAYF6UbmkO/iBOUOLbD7FYEX14zxe86cHM9PPzeHxvZufvfeEZ/JN15wVSEcCsRhtwy8RWv9sYtjLxjJOX1NV28fJ+vauHJB5qDtX1mdx6t7ynn+01Pcec74q53r6O7jsU0SxzBWthYDHKlqYdX0NNbvqxhUVerIgEUwfHxifFQ4C3IS2XyslnsunulXWYOVI1UtpMZGkGpiUS2emszT20rYWdLARXMyPD6nY+bQijzfDAkqa+xgWqoomAy7i8qbzKHCyhbS4uT3TI2L5OYzp/L3rSe5cflUZk6Od38CE3aebOCHL3xGXGQYqbERpMRGkhYXQU5yNDefOc1tFbY/cKUIKrTW/z5mkgSA4tp2bJpBFgHAgpwklucm8+SWYr60MrffbTReeGZ7CbWt3aw5I51NR2ro6u0jMswz83ui09TRQ0VTJ/OyEtlZ3OBUEbiyCABWFqTy5w+LaO3qJS5yQnVU8QlH7D2GzFg8NQmA7l6bRzUEBtlJ0USFh/gsYKy1pqyhg5UFaYCDReBFvyEJiA/8nvdcNJM39lbws9cP8PSdZ5qO33THXzafoKa5iynTYqhp7aKwsoXatm66e21MSYnhvJljn0jpNkYwkenPGEofbt5+ZXU+pQ0dvHtwfGWPdPb08ecPj3N2fiqX2ytha1utnHeDow4ui6ykKKeuoZqWLiLDQpze5FcWpNFn03x6ot5vsgYrWmuOVbeauoUA8lJjSYwW15ynGUMgYyvz03yXOdTc0Utbd1+/ayg2Moz4yDCPLQKbTXOkqnXQyj85NoJ7Lp7JluN1Iyo+bWrvYcPBKtYtzeFvX17BG98+hy33X8Cef7uIsBDF9gB931wpgnHhvvEnxhfOTBFcNCeDaakx/M84G2/3/KenqGnp4tsXTO9fzZ5O7qFdJQ1UNDmvDzDaHpyREU9WUrSLGIFUFTtb0S2dlkxEWAhbjteOXugJRnlTJ61dvf0VxUMJCVH9VoE3MQKQNG5fFZWZpa9mJEZ5HCw+1dBOR09ffwsNg5tXTGV2ZgK/XH/I695Ir+8tp7vXxrolOYO2x0SEMS87cfwpAq31hF8KHatuJTspmuiI4W6V0BDFl1flsbukkZ0nnf8pKps6+fo/dlI9BkHbrl6xBpbnJnN2fupppwia2nu46bFP+I83Dznd50hVC7ERoWQnRZOdHE15U6dp7rdZVbEjUeGhLJ2azOZjgQ0YN7Z3c98Le9lf1hRQORwZCBQ795Eb9QTeWAQgPYd8NbbSMXXUwJtaAiMN+YwhsYDQEMXd5xdQ2tDhdULBi7tKmZkRz7zs4dN+z8xLYW9pU0DSlseX83uMOVbdOiw+4Mh1y6TA7PGPnJeX//zNg7x9oJItY5Bh8sLOUiqaOvnOBTNQSp12iuDl3aV09drYcqzWaWFPYWUL0zPiCQlRZCdF091rMx1GYlZVPJSVBakcrGgOWLuJ2tYubnzsE57fcYp/bD0ZEBnMOOoiY8jgumU5fP3cgmHppe4omBSL1lBUO3qrwGxMZkaC520mXCm8C2dnkBAVxou7Sj2W53hNK7tLGlm3NNvUEl2Rl0J3n409pxpNjvYvp60isNk0RbWuFUFMRBg3nzmVdw5UUlI3vMBsy/Fa3txbAeD3yVbdvTb+tPE4i6cmsXq6BL9SY08fRaC15tntpwgLUTS093Co0rxI/Wh1CzPtN6isRLkBlDcOv/DdWQQAK6enAvBJ0dhbBVXNndzw6FaK69rIT4/l0+LxY6AfqWolPT6SpJjhGVcGWUnR/OiyWV5nwAyMrRz99VTW2EFEWAipsQNyTk6MpLqlyyOXTmFVKznJ0aZxpKjwUK5amMVb+ys8TjN/cWcpoSGKqxdlm76/bFoKShEQ99BpqwjKGjvo7LG5VAQAX1qZS2iIGtbArKfPxs9eO0hOcjQZCZGcrPOvInh5dylljR18Z+2M/tVERFgIyTHh1LRO/FqCXSWNFFa1cPf50wFMTfLa1i5qW7v7V3CGS2BonKC3z0Z9e7dbi2BBThKxEaFjYu05UtrQzvWPbqWyqZMn71jBDcumUFTbNmwwfKA4atJawlfkpcWiFD6ZVlZmryEIcVBGkxOi+tuLuKOwsnlYfMCRdUtz6Oyx8dY+90HjPpvm5d1lrJmRxiQn09oSY8KZNTnBUgRjiRGQcqcIMhKiuGphFv/ccYqm9hkPcsUAACAASURBVAHN/4+tJymsauEnV8xh+qQ4ik0sBl/yxMcnmJ+dyHkzB7fbmBQfRXXzxLcInt1eQmxEKHetySc/LZbNx4YHcR2Lf2DAPz1UEdS3daM1pJvUEDgSHhrCirwUNo9hwLi4to0bHv2EhrZu/nHnmZyVn9qfU7+jOPCVzkYmjbcuH0+JCg9lSnJM/4yQ0VDW0EFW0uCbboaHRWXdvTaKatqGxQccWTwlifz0WF7Y6d49tOV4LRVNnaxbmuNyvxW5yew82UDPGM9FsRSBScbQUO5cnU97dx/PbJcJZrWtXTz03hHOmZHGJXMzyE2N9atFUNfaxZGqVi6fnznMt5geH0mNB6ubsaSls4f3DlbxwGsHuOHRraPu29TU0cMbe8v5/OJsYiPDWDk9le0n6oddLEcdMoYAEqLCiIsM6++CaVDjpobAkZUFaRTVtI3J1LLWrl5ueGwr7d29PPPVs/oHvczLTiQ6PDRgGSWOlDV20NHT5zJQPFoK0mN9YhE4VhUbTE70rM3Eido2em3apUWglGLdkhy2F9e7vf5f3FlKQlQYF852XWC3Ii+Vjp6+MU8OOG0VwfGaVlJjI0iOdb0qBJiTlcCq6ak8ueUE3b02Hnz7MB3dffz0qrkopchNjaWhvWeQxeBLjPa3y3KHD3NLj48cFzGCzp4+/vj+Ub7wyBYW//sG7vz7Dp7ZXsK2E/V8bLJ694bX9pTR2WPj5hVTAVhVkEZbdx+fDQmqFVa1kBgdziT7DV4pZVpL4KrP0FDOLpA4wdYi/1sFW4/XUdXcxUM3LGJedmL/9vDQEJZMSxoXimDA6vKPawjESj8xyrGVXb19VLd0DcoYgoGiMncB48Ih1qUzrl2SjVLw4i7nQ5BaOnt4+0AlVy3Mctt3aXmeXONjHRM6bRXBsepWCty4hRy585x8qpq7+OX6Q/xzRylfXp3X71aaZp+HerLeP1bBzpIGwkMV8x1uDgaGIvBlx8aR8PHRWn797hHau/u4a00+z3z1TPb+9GIiQkMoHkUgXWvN09tKmJed0H9zPLsgFaUYltp5pLKFmUP642cnRVM+pO7AKMDzxCKYk5lAUkz4mKSRbjleS2RYSL/ycWR5bgqHKptpDnD/qx32Rcl0P7mGQALGXb02pw0DPcFY8Q+1CFLjIgkNUW5dQ0cqWwgNUeTb+x85IzMxmtXT03hpV6nTTLb1+yro7LHxBTduIRBXb35a7Jgr/dNSEWitOVbTalpI5ozzzkhnxqQ4ntxSTHp8JN9eO73/vdw0+bL4K06ws7iBedmJpquJ9LhIunpttHR53lrXHxg523/78gruvXQWKwvSiAoPZVpqDEWjUASflTZxuLKFm+zWAMgQk7lZCYN891prCqtamDEkiClFZYMvem8sgpAQxdn5qWw9Xud3Zbv1eB3Lc1NM24WsyE1Ba/kuBIq/bSnmkQ+Oc+HsjP7KYX9gLNCOVo+8xXqZSeooSA3ApPhIt20mDle2kJcW61Hrli8szaG0oYNtTm7eL+4sIz89lkVTkjySfUVeCttP1I/p7IPTUhHUtXXT2N7jNlDsiFKKr9ob0N1/2SziowYuhKkpdotgBDc8dzeXrt4+9pY1sXTqcLcQMG5qCcqbOogIHZyqB6IkR2MRPLuthJiIUD63MGvQ9lUFaewuaejvLW/0xx9qymclRVPf1j2oB31taxcxEaHEethDaGVBKmWNHZT4cUZFbWsXhytbTK0BkEZuYSGK7QFKI3144zF++toBLpmbwcO3LPbrZ83NSrBXdY/cCjOsCbOCNk9qCY5UtbiMDzhy8ZzJxEWa1xScrGtje3E965bkeNyXaEVeCs2dvf3uqbHAb4pAKfUXpVS1Umq/w7YHlFJlQ9pSjzmeZgwN5bplObz2rVVcs3hwHnBUeCiZiVFeWwRtXb2s/K/3+cfWYqf77C9rprvXZhofgPGjCCqbOslIjByUqgeSDniyvn1Eq5uWzh5e+6yczy3MGqR4AVZOT6OnT/OpfYVc6KT4JztpeC2BJzUEjpxtb1o2mhuTu4C5kQ670okiiI4IZX7O2Lcg0Frzq7cP8//eKeSaxdk8fPMSvzc4jIkI46z8VDYeNp1Z5RHG/3uyyVCcyW4UgTHq1NPuotERoVy5IJP1+ypos1vmfTbNP3ec4vpHtxIeqrh2iXntgBlGlthY/q/9aRE8CVxqsv0hrfUi+2O9Hz/fKSNVBEopFuQkmWr2aakxXmcOGZXCL7gINO2y+2SXTBvfiqCisZPMxOGrr9zUWLp7bcP89J7w6p5yOnr6BrmFDJbnJhMeqthiD0QfqXSiCExSSD2pKnakID2WjIRI05RVT/jwSA3nPLiRj486P37L8TriI8NM40AGK/JS2Fva6LQFwTsHKn06sc5m0/zbqwd45IPj3HLmVH5z3cIx68S7dmY6RbVtI7YmyxrbmRQfaaq0Jie6bjMxNPvME9YtzaG9u4+391fy8dFarvzDx9z7wl4yE6N57q6zTa8NZ+Qkx5CVGDWmisBv/XW11pvsQ+vHHceqW4mJCCXLwxF6npCbGst7hzxfwdhsmr9uPoFS8NmpRiqbOk1XLztO1jM1JYZJ8eaypseND0VQ3tTBMhNllZsmbrPi2nZykt0PMTfQWvPcpyXMzkxgQc7wm2NMRBiLpyb3xwmOVLWQFhdJyhDXlFlRWU1Ll9sgoCNKKVYWpLHpSA1aa69bDxt55v/ccYrVM9JM99l6vJYz81Nc3mhX5Kbw6IdF7C5pHOZCKqpp5RtP7eSqhVn8/kbvXDdaa/6yuZj3DlbR3i0dO9u7emnp6qWls5evrcnnR5fNGlHL5ZGydlYGD7x+kPcPV/Pl1XleH1/e2DksY8ggI2FgZGVMxPBboNFjyJt5A8umJTMtNYZ/e3U/bd195CRH84ebFnPlguEp354g9St1I/q+jYRAxAi+pZTaa3cdmS9z/cxxe6DYl3/gaamx1LZ20eph0HZjYTXFde18Z+0MADYcHF6dqLVm58lGljqxBgASo8MJD1UBrSWw2TRVzZ1MNln15NkD6Se8tJb2ljaxv6yZm1dMcfp/WlWQxoHyZhrbu8Wna5LSmBEfSYhi0FwCby0CgLPyU6hr6/a69UFrVy8bDlYSHqp492Cl6fejrLGD4rr2fheUM4wWBGaphb//11FsGk6OwD159zO7+PkbB2nu7CE5NoKZGfGsmp7GuiU5/Oa6hWOuBACmpsZQkB7LxsKRuYfKGjucNrzLSJD/vbNagsKqFqLCQ/pjf56glOIO++yS+y+bxXs/OJerFmaN+O+2Ii+VmpYuvxeqGoy1IngEKAAWARXAb5ztqJS6Sym1Qym1o6amxqdCHHfTbG4k5BoppB7e8J74+ASZiVF8a+108tNiTecelNS3U9va5VIRhIQo0uICW0tQ29ZFT58eVsUJkBEfRVR4CCe8vIE+ve0kMRGhXL3YuW911fRUtBb/urNBKWGhIUxOiOpXBD19Nhrae7yKEQAsnSZ+W1edaM14Z38lnT027r1klr0dQcWwfQz31qrp5vEBg8SYcGZmxA9TBEeqWnjts3IiQkO8Cmifqm9n3SNbeHt/JT++fDZvfHs1T96xgodvWcL/u24hD3xuLuuWeh7k9DVrZ01iW1F9v9/dU7TW/e0lzHA3svJIVQszJsV73Sfp9pW5fPbTi/nauQUez2l2xkCcYGzam4ypItBaV2mt+7TWNuB/gBUu9n1Ma71Ma70sPT3d2W5e09bVS3lTp88VgTEOz5MV2aGKZrYcr+P2lbmEh4Zw8dzJbD1eN6wgzSgkc6UIIPBFZRX2wJyZHzQkRAruir2wCJo6JEj8+UXZw4LEjiycIr2Ant9xymW1a3bywFyCOi9qCBwpSI8lOSa8PzjtKa/sKWNKSjRfWZ3HtNQYXt49PB609XgdqbERnOFBbv6ZeSnsPNlAr0NV9e/eO0JsRBh3rM6lvq3boyZom4/VctUfP6bC3s/oq2vyA3bDd8b5sybR3WfzuiCxtlWmfTlTBO5GVhZWtoyoctqXf7+C9FhSYyOcpqT6mjFVBEopx+HA1wD7ne3rLwaG0XjuI/YEo6jMkxveXzefIDo8lBuXTwHg4rkZ9Nr0MDN4x8kG4iPD3H4p0wNsEVQ0GYrAPI6R52UK6cu7SunssXHLmcODxI4YvYA+KBSL0dnfKSsput8i8KaGwBGlFEunpfQrZ0+obulk87Farl6UTYi96+TWorpBg3W01mw5XsdZBanDMq7MWJ6XQnt3HwfKpfvqgfIm1u+r5MurclmYI3nq7qyCF3eWctsT28iIj+K1b61izRm+W2j5kuW5KcRHhnmdPVRuMofAEVcjKxvauqlu6fJr5bQnKKX66wnGAn+mjz4LbAVmKqVKlVJfAR5USu1TSu0Fzge+76/Pd8ZIM4bcERsZRnp8JCdrXV+Eta1dvLKnnHVLs/vb+C7KSWJSfCTvDokT7DrZwKKpSW5N1ED3GzJubM4UQW5aLCX17YNWsc4wKokXTkka1GbBGaumD/jVnXXEzEqKprKpkz6bdjur2BXLc5M5UdvWfw53vP5ZBTYNn7e3Hb5mcTZawyu7y/v3OVHbRmVzp9O00aGsyB2cWvjQhqMkRIXxlXPy+33a7lJVH//4BLMzE3jpmyv7LdnxSHhoCOeckcb7h6u9KuYbGEhj/n10NbLSWRpyIFiRl0JpQ4fTKXu+xG+KQGt9k9Y6U2sdrrXO0Vo/obW+TWs9X2u9QGv9Oa31cIepH9Fa89b+SiLDQvxyAeSmxri1CJ7ZVkJ3r40vrRzIhAgJUVw0J4MPCmv6UwObOnoorGphmd037YpJ8ZHUtXrWY90fVDR1EhEWMixjxyAvNZZem3Y6SN6RT4sbOFrd6tYaMDAGk2clRjl1I2UnRdPTJ0rAsAjSvbQIYKDXk6ddQF/dU8a87IT+RUduWixLpibx8u7S/hubUZuwyk2g2GBSQhS5qTFsL67ns1ONvHeoiq+ek09idDhT++NUzhWBzaYpqmnl7PxUjwvqAsn5MydR3dLVbwF5gnHjzElyHux1NrLSSL+dNXn4BLGxZrld6Xvzu4+U06qy+Nntp/j/2zvzMKmqM+H/3urqrXpfWRp6AaSbJiACEhcUQUDjMokxGv3yKdk+k895MonLTDLjmJnM6mgWzTpxEkdmRs3ixGgUNSyZuCECAiKCItA0Te8rvVWvZ/44t4rq7qruqqarblX1+T1PP1V1763q89a9dd9z3nXrew3cvWEhiWGIhy7JSxv3R9g3OMR/vnmSK8oLxqxIrlo8k57+IW+s+v5T7Sg1sX8A9Ox2WOnyynZQ1+FmVlZKQBuppwRHMM17nth1kowUJ9cvnT3hsQAVMzPIS0saN9TPYys+3d7rXTmFahoCXQU0yekIymF8rKmLd2o6xjQhuWH5HD5o6PL+uN841szsrBSvaTEYLizNZXdVK9/Z+gE5rkQ+Z4VXZqYkkuNKHNc0dLq9l77B4ZDqbNnJFeWFAOwIwTxU09ZLWlICmamBFV2glpXv13eSmeL0RhbZyaJZmey7fwMbKsevWDoVTBtFcLShk797/hCXnZfvLRUx1ZTmuag/4w7Yb/WFd+po6uzj85eOjYu+aF4eGSlOXj6kzUN7q1pxCCwrnrg+icfMYVfjkrr23oBmIfDNJRhfEbR09fHiwXpuXD7Hbx9pfzgcwk9vW8F911YGPMZjKz7d1ktTZx/pyc6gP9+XZGcC58/JCsph/Oy+0ziEMaUxrlsyi8QE4Zl9pxkeVuw81sLF8/NDcjSuKsulvWeAVz5o4ktr5o/ooFWc6xpXEZz1kcWGIijISOb8OVkhKYJaK3R0vO/UX5mJTvcA2w83smROVlQ4zhMcElR15KlgWigC98AQX3lqH64kJ9+56fygnHKTwRs55KcKqVKKn792ggWF6VzmJ6koyelgXUUh2w43MjSs2FvdRsXMTL9t8kZjd3ZxXYfb2xbSHwXpyaQlJUy4Inh6bw39QxM7iUezsjR3XJ+Px1Zc2947qRwCX1aU5HKodvwG40opfru/lkvmj+1GlZOWxNryQp7dX8uh2jO09QwE7R/w4AktzE9P5vaLS0bsK85LG1cRHLfCeKc6WCKcrK0o5EBNe1BdxUCvegI5ij34a1n5wItHaOx08+dXVZzTeGORaaEIHnjxCEfqO/n2TUsDtombCkotRVDlx2F8qPYMh2rPsOmS0oCzjY2VM2nt7mfXiRb2VbcHrC80moJ0LZMdimDISiabFcAxBzoCoqwgjRMT2K6ffKuaVWW5nDfFjrqMlEQyU5zUtvdadYYmP8u6sDSHgSE1pheCL/tOtVPd2sPHl/k3b31yeRHNXX08+PIR4Gxv5GApznVx9eKZ3H/dojGZscW5qZxu6w3omD/W1EVWamJAf040sq6iEKXwRoeNx1snWjne1B0wdNTD6JaVu4638MSuaj5/aVnQVULjibhXBDuONPD4G1V89pJS1lWE19ZWPE5S2XMHanE6hOuXzhqzz8Oa8gKSnA6+v/0oPf1DQfkHAPIz9I/ajsih5q4+BoeV36xiX0rzxg8hff1YMydbekJeDQSLDiF103TOKwLLYTxOGOlv950m2eng6o/M9Lt/bUUhWamJvHq0mXn5aSHVoQGtWP/1thXeaCRfSnK1Y74uQNaszqpPiwrTR7B8ZHYW+enJ7Bgny7inf1B3xHt0J/kZSWy6pHTcz/RtWekeGOIbvznI3NxU7t64cCqHHjPEtSJoPOPm3l+/w6JZmXzjY+Ff7nlmWqPTwoeHFc8fqOXyhQXekFF/pCc7Wb0gnzePa2dksIrAlaRbMtqxIvDGbE9Qt6ksP42ath76B/3PVJ94s5rctKSAN89zpcjKJWgOsfLoaLJdSSwoTGdPgHLQA0PDPP9OHesrZwSMYkp26mqVQMCy05Nlbu74kUPHmrpjxj/gweEQ1pYX8MoHTX57+b55vIWrH36Vx9+oYtPFpbz01csnDP/0bVn5yPajnGju5oFPLvVbe2g6ENdS/9OWw/T0D/KDW5edc8p3sPirQrq3uo3aDjd/cfXEymhj5Qx2HGlkRmbyhMtbX+zKLq7vCJxV7EtpXhrDCk619Yy5ETV39bH1cANfXF0WthLHRTmpvHm8he7+oUmFjvpyYWkOL7xTx/CwGuNv2vZeA63d/WOihUZz08q5PLGr2hsVM1V4VqX+/ARn3AM0dfbFTMSQL+sqCvn13hru+dUBsl1nFWxLVz8vHKyjJM/FL++4iI/OC06xepLKdhxp5Nd7a7h55ZwROSnTjbhWBPdfV8knLigKa1u90ZTmjW0z97sDtSQ7HawPIgxsfeUM5JmDVoGx4Jfv42UX//SPx3jx3Xru3rBwyrNIayfIKvbg7eLWPHZG+vKheoaG1bh1hc6V2dmpdFvRXPnnYBoC7TB+6q1THG3sGhG2OjA0zEMvv8/8gjTWlo//PS+bm81rX18bkrIPhpmZKQFrDnkcxfPyY8dR7OGyhQUsnJHOq0dH+gkSHMIXVpdx78bykCLBPC0rf7H7FAUZydx3TeCos+lAXCuCvPTkKZ9xTURJnovf7j+Ne2CIlMQEBoeG2XKwjvWLZgQVAZSfnsyDNy5l8eyJs2p9KchI5nC9/8STp/fWcLSxi9sfe4u15QXcd23llGVW17X3kpLoGDFL80fZOLkEWw7WMS8/jYoQyv6Gim8UyVSsCECXCPdVBE/uquZ4czc/37QyqLr9oZTlDpYEhzAnJ5VqP5Frx6ys+lhcEaQnO/n9XWum7PM8LSvrOtz8/ccXkzXB9RvvxLWPwA5K89JQCmra9Ixs5/EWmrv6uf78wE7i0dy0ci6Vs0PLbAxkGmrr7udoYxd/duV5/OXHKthT1cZVD7/C3z53iLYpSECrO6NDRydaveS4dOTOaEXQ0tXHzmMtXLNkcnXbg6XIJ6rpXFcExbku8tOTR2QYn3EP8PC2D7h4Xh7rKiI7+RhNcZ7/XIJjTV04HRJSeeV45uJ5eXxqxRyu/kjwv814Ja5XBHbgLT7X3MOCwgye219LerIz7CuTgoxkOt2D3pWIB090y6Xz8/jovDxuXDGH7239gP/YWcWrR5vYdveac7oB17X3+m2oMxodQpo+pgTHS4fqGVZwzZLw/hiLfMoNnEvUEGhZLizNYY9PhvGP/vAh7b0D3HftItsjcopzXew92Tamqcnxpm5K8lxhyaqPRb776WV2DyFqMFfEFOPNJWjppm9wiJcO1bNx8YywO6sDdSrbU9VKUoKD863Y6Pz0ZP7xhiXcf10lx5q6OdV6bgWtdHmJ4OzcZXmuMTkWWw7WUZafxqJZ4fXjFGQk47Qcu3lTEEO/oiSHU629NJxxc6q1h39/vYobLigKqlBeuCnOddHpHqSjd2Q56mNNXcyLsYghQ2QwimCKybZMICdbevjj+010ugfHlBkIB97s4lG5BG9VtbJkTtYYRXSRFV2xt3ryZW4Hh4Zp7OwLWOVxNKX5adR29Hqzcs+ahWaGfRad4BBmZqWQmeKcEqW80ioItqeqjYdefh8B7t1Yfs6fOxUU+wkhHRwapqol9kJHDZEhnGWoHxORRhF512dbrohsFZGj1qMtrSrDiYjo+vst3fzunTpyXIkRCUvzV2ait3+IgzUd3iqGviyckUF6sjOk+vqjabIqngZjGgLtMFbqbGjjy4caImIW8jA7O/Wc/QMeFs/OJCXRweY3qnjuQC3/77J5E5Y1iBSeUie+foKatl4GhlRMlZYwRI5wrggeB64ete0bwHal1HnAdut13FGSl8b79Z1se6+Ba5bMiohNttCPIth/qp3BYcWqsrH6NsEhLJubzd6TgUslTESt1ZlsvDpDvnjMZh6H8ZaDdZTmuaicFZmSv3deMZ+71k9N5mhigoNlc7N5q6qV/PQkvnzF/Cn53Klgbq4+H76KwFtsLgYjhgzhJ5z9CF4BRtsdPg5stp5vBj4Rrv9vJ6V5Lho7++gdGIqIWQggNy0JkZGKYHdVKyKwoth/T4PlJTm8X3/Gb0N1D+6BIb/N0sEnmSwE0xBoRdDa3c/O4+GPFvLlivJCrp/C8+FZad21YWFQocGRwpWkmyRVt/hRBPlGERjGEmkfwQxPMxrr0d44uzDhWZrPzEzxa5YJB84EB3lpSSN8BLurWimfkREwRnpFSQ7DinELqP3bK8e56V93+o3/P9uZLLgVQVZqInlpSVQ1d3uTyCJlFgoHt6wq5p4NC/n0yrl2D2UMxbmuEVVwjzV2k5+eNO3j5Q3+iVpnsYjcISJ7RGRPU9PEVQejCU/9/euWzgpbyWt/5PtkFw8ODfP2ybZxK5gum5uNCOP6CbYebgBgu/XoS227G1dSApkpwc+GS/PTONHczRarLMDiEPMloomi7FS+cuV5QSWPRZriXNeIiLDjzSZiyBCYSF/BDZ4G9tZjwHKCSqlHlVIrlVIrCwqis7l2IBbPzuKWC+dOWAFxqinI0DXWAQ7XddLdPzTuiiQrNZGFhRkBFUHDGTfv1HQA/jtE1Z/pHbczmT9K89I4XHeGNyKQRDadKc51UdvRS9+gjtDSxeaMo9jgn0grgueATdbzTcCzEf7/ESElMYEHblzqrQQZKQoykmm2FIHHru9pYhKI5SU5vF3dxrCffsfbD+ub/7qKQt460coZ98i49Np2d8iRMmX5Ls64BxkaVlwbw2ahaKc414VSuitba3c/rd39JnTUEJBwho8+BewEykWkRkS+ADwAbBCRo8AG67VhivCUmVBKsbuqlaLs1Ant9ytKcuh0D/Kh5Uz0ZfvhBubkpPLlNfMZHFa8drR5xP66jl5vFcdg8TiMi3Nj2ywU7Xgy3E+29nA8xtpTGiJP2EIdlFK3Bth1Zbj+53SnID2Z/qFhOnoH2F3V5rcl5mg8PQ/2nmwbUcO9t3+I1z5s5tZVxSwvziYrNZHthxu9zt0BK5lsVsgrAq0Irl1qzELhxJNUdqq1hxSrtLdRBIZARJ+XyzBpPEllu6vaaO7qCypiqTTPRW5a0hg/wWsfNtM3OMz6RTNwJji4oryA/3m/0WtCauzsQ6mJG9KMZtHMTP7i6nI+f2lZSO8zhEZBRjIpiQ5OtvRwrKmLJKeDopzoSHgzRB9GEcQRHkWw5WAdcLZc8niICMuLc3h7lCLYfriBjGSn18ewrqKQlu5+DtToUNM6qzNZsFnFHhwO4c4rFpxz4TfD+IjoKqPVrVoRlOWlkRDBCDZDbGEUQRzhyS7e9l4DOa7EoHsOrCjJ4biV5AW6tea2w41cbvVQBlizsACHnI0e8jSkiZayCoaxFOemUd3Sw/GmbuYXmoghQ2CMIogjCtL17Lyzb5AVIXQ48/gJPKuCd0530NzVx/pFZ/P9sl1JrCzJ9SqCem8yWWgrAkPk8CSVnWwd2x7UYPDFKII4IjPVSZKV3OSvvlAgls7JwukQ3q7WimD74QYSHMLaUT0U1lYUcqj2DPUdbmrb3aQnOwM2aDfYT3FuKu6BYYaGFfNMDoFhHIwiiCNExGt7D6W0RUpiAouLsrwO463vNbCiJIds18i6/VdaK4QdRxqp6+g1q4Eox1PqBEzEkGF8jCKIMzzRIqH2PF5enM2Bmnaqmrs5Ut85wizk4bzCdObkpLLjSCP1He6QQ0cNkcU3odGUlzCMh1EEccZl5+VzwwVFXidvsKwoycE9MMwPdnwIwPpFM8YcIyKsqyjk9Q+bqW7tCTl01BBZ5uSkIqKLH0ZTdVRD9GGujjjjnkl2yfI4jJ/ZV8O8/LSAM8h1FYX8x86T9A4MhRw6aogsKYkJzMxM8SbxGQyBMIrAAOhS0rOzUqjtcHt9Af64aF4eqYkJ9A4MBd2QxmAf37yukpwp6NFsiG+MacjgZbm1KvBnFvKQkpjgbb0ZbEMag318bMksb39qgyEQRhEYvHxiWRGrF+R7zUSB2LhYK4qSXGNyMBjiAWMaMnhZDJEQyQAACo9JREFUXzmD9ZWBVwMePrV8Dh+ZnUVxXmTLbBsMhvBgiyIQkSqgExgCBpVSK+0Yh2FyOBxCpSkhbTDEDXauCNYqpZonPsxgMBgM4cT4CAwGg2GaY5ciUMDvRWSviNxh0xgMBoPBgH2moUuVUrUiUghsFZEjSqlXfA+wFMQdAMXFxXaM0WAwGKYFtqwIlFK11mMj8Aywys8xjyqlViqlVhYUFER6iAaDwTBtiLgiEJE0EcnwPAc2Au9GehwGg8Fg0IhSKrL/UGQeehUA2jT1pFLqHyd4TxNwcgr+fT4QD5FK8SIHGFmilXiRJV7kgMnJUqKUmtCkEnFFYCcisicechbiRQ4wskQr8SJLvMgB4ZXFhI8aDAbDNMcoAoPBYJjmTDdF8KjdA5gi4kUOMLJEK/EiS7zIAWGUZVr5CAwGg8Ewlum2IjAYDAbDKIwiMBgMhmlOXCoCERG7xzBVxIss8SIHGFmilXiQRUQSrMeIyhI3ikBEFovIFQAqxh0f8SJLvMgBICLlIrIE4kKWuDgvIrJaRH4iIndCzMtyqYhsBv5aRHIjLUvMO4tFxAH8EFgHVAO7gGeVUntExKGUGrZ1gCEQL7LEixwAIuIEfgqsBuqA3wG/UkqdEhGJpZtPnJ2X5cBm4BHgE8BRYLNSar+tA5sEPtUWvgdcDvQCW5RSL0RqDPGwIsgBMoBFwGeAFuAeEUmPpQvbIl5kySI+5AAoATKUUuXA/wcKgDtFJDWWlIBFNpBOfJyXVcBupdTPgC8CPcA1IpJv77AmxQrgsFLqceAeYD9wnYjMjdQAYlIRiMgGEdlgvcwELgZcSqkm4L+BVuBPrWOj2m4oIp/yLG2JYVlE5JMi8j3rZR4xKgfo2aaILLReJgIrRSRRKXUYeA5IA260bYAhICJlIpJivcwFLiEGz4uI3Cwid4vIJdamt4F0EZmplKoHdqBr8Vxq2yCDREQu8rm+AHYDc0RkrlKqDXgdaAduiNSYYkoRWLbNXwB/BbQBKKVOoL+4r1mH1QG/AS4QkdnROmsTkXQR+W/gXqBNRJyxKIuIVIrIk8D9wFetcX4I7CSG5ADvTfMF4EfAf4rIBqXUEWA78H+tww4A+4DzRSTbpqFOiIiUisiLwM+AJ0Sk0jovrwB3W4dF/XkRkQQR+SbwdWvTT0XkeqAbqALWWNv/CHQAc633RZ1SE5Fs6/raCtwsIunWLjfwGnCz9fp94D0gz0eJh5WoVwSeEyoiueiLuFUptVYptcfnsMeBS0WkTCk1CDSgv9zUSI93PEZdnHOBBqXURUqpp4Aha/vjaFnmRassPufkcuDfgDeVUheg7bUftQ77ObF3Tu4F9iulLgaeBW63tr8KXGzdLLuBGmAO2pYbNfiRZZdS6krgD8C3RKQSfX1dFM3Xly9KqSGgHLhHKfVd4FvAV9CVi+uAZZaSG0TfQG+w3hd1Sg29knwZPf40tD8AoAl4E1giIqssmU+jG3i5IzGwqFcEQAqAUqoVeAhIBhCRz4rIVSJSopT6A3qW9pB17Lto226fPUMOiK92X4q+mWCZhv5GRFajZwJvAN+GqJXFc9N4D9iolPq+iCQBCwCPrfkAevn+IEStHGCdE+sm2g0MWNszgaMiUoKegDQCf27t2w4UWcdEEx5ZPJ0HDwEopX6ItqnfAtSiTRFRe15E5HYRWeOz4moAcqxV89PAMWA9+jy4gX+wjisCdvvIbzs+smQqpU6jy0T8Cj3uVSJSZN3430Tfw75nrRQWA9Ui4orEOKNWEVh+gK3AQyJyi7X5EeBCEakD/gS4BvidiMxHzxSKROQHIvIuun9BRzQsEX1keVBEbrU2vw3UichjaHt6O3Af8HHgu0ChiPwwmmQZJcctSqlmpVS3iKQopfqBg2gnJJat8+/Qts9oPicPicjN1gzyNeA8EdkHXI2edf4SqECvcNZbfpCDaEXXac/oR+JHlkG07f8CETlfRM5HN38qAxLQN86oOi+imSUifwA2oa+jH1k3xWZgCdrRDfo+cBvQqJT6FtBumVxuAX5myW8bAWT5iYjkK6XcSqkeYBs6OGQdgFKqXin1CDqS6zG0KfJfrGPDj1Iq6v7QM8td6JviBcATwF9Z+64HNvkc+5j1hQHMQDvD/sRuGSaQ5R70TeY7wF4g0Tr2NuBR63lhNMniR47/8jknnvGvsbYX+LyvIJrkCCDLk8C91r5y4Dc+x34T+L71vNS6/j5ptwzjyPIUcCc6aut+4Hm0gltpyfk1631R81sBEqzHhcB/Wc+dwI/RCjgbbVK5HO3oBj2rvstz/flec1Eqyw98rytr+11opZyFjkwDragzIj5uu784ny/FATis558Bfuyz7/PoGXOh7/HW443AT+wefwiyfMGSJdu6sHcA/8fatxT4ree9dv9N4pysR8fZO+0e+yRlmYFWXI8Ai6x9q4Gno+WchHB9FViv5/ns+1Pgi9ZziQI5nMA/Af+CnkRcj84F8JWzEW0m2WQphk9b+54APmq3DCHIImifxhqfbenAw8BbaPPXbLvGHxWmIRH5HNoB9/fWpoPArSJSar1ORNsFv+15j1JqWEQ2AX8DvBSxwU5AELI4gRPAg0qpV9AXwj0i8nXgF+jZm+1RD5M8J9vQM89LiCKClOW4tb8THWb5ZyLyVXQy2TYgKpyPQV5fx9DJSaCvNUTkDrSSeBvsd6aKyBr0ajgH+BAtzwCwVkRWgf6No02+DymlNgO/B263THdOtOy2E6QsCm0q/Vuft16LXr0dAJYopWojOOyRRIEmTUfPgr+KvkgrrO0Po5e5r6PNDUuAF9Amkzy0Y/h/gAvtlmGSsmwBZlr7LwS+BFxstwyTPCceORKBO4BSu2WYpCwvoqM5FqEjOzYDF9ktwzmclxnW/q+hHcTR9Fu5DLjN5/WP0Ql7nwX2WtscwEz0imyutW0mPqucaPgLUZZfeX4faHPe5XaPX6koMQ0BxdbjA8AvrecJ6JnZauv1XHTom9P6K7F73FMgS4rd450COf4dSLZ7vFMky2Ygye7xTuH1lWy9dtk9bj9yuNARgB6b+meAf7ae7we+Yj1fCTxl93jjXZaoMA0ppaqtpw8DZSJyldIhVR1KqdesfV9Gp5GjlBpUSp20YagTEqIsA/4+IxoIQY5ewNYojYkIQZZuzuZzRCUhXl+D1nsiE3kSAkqpHqVUnzV2gA3oeHqAzwGLROR59ErnbTvGGCyTkcVu0+8Y7NZEfrTrl4A/+rxehU7u8ZpSYuUvXmSJFzmMLNH3h17NONBmuQXWtgXoYIrVQJHdY5wOskRV9VGxKiCKyNNoD3sf2lF3VCl1zN7RhUa8yBIvcoCRJRqxZsZJ6FIYz6AjuFrQ5pQzdo4tVGJZlqgwDXmwLmwX2iF8K1CtlHopli5sD/EiS7zIAUaWaETpmegFaLv63cAzSqlN0X7j9EcsyxI1qdg+3Im2o21QSkVN2vskiRdZ4kUOMLJEIzXorPrvxrgcEKOyRJVpCM4uee0ex1QQL7LEixxgZDEY/BF1isBgMBgMkSWqfAQGg8FgiDxGERgMBsM0xygCg8FgmOYYRWAwGAzTHKMIDAaDYZpjFIHBYDBMc4wiMBgMhmnO/wJUIgSFtrWcewAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2aaacffcba20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "out=figdir+\"VIKING20-K301_NA_transport.pdf\"\n",
    "\n",
    "nac.plot(label='NAC')\n",
    "gs.plot(label='GS')\n",
    "plt.ylabel('Transport (Sv)')\n",
    "plt.xlabel('')\n",
    "plt.legend()\n",
    "plt.savefig(out,facecolor='w',format='pdf',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create dataset\n",
    "out=dir+\"VIKING20-K301_1y_1949_2009_NA-transport.nc\"\n",
    "\n",
    "nac.name=\"nac_transport\"\n",
    "gs.name=\"gs_transport\"\n",
    "nac.attrs['units']=\"Sv\"\n",
    "gs.attrs['units']=\"Sv\"\n",
    "\n",
    "data=xr.merge([nac,gs])\n",
    "data.to_netcdf(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:        (time_counter: 61)\n",
       "Coordinates:\n",
       "  * time_counter   (time_counter) datetime64[ns] 1949-07-01T12:00:00 ...\n",
       "Data variables:\n",
       "    nac_transport  (time_counter) float64 12.44 14.04 21.27 18.66 16.08 ...\n",
       "    gs_transport   (time_counter) float64 27.85 30.02 27.39 29.64 34.09 ..."
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SSH shkpwagn@hdata2.hlrn.de python_py3_std_/gfs2/work/shkpwagn/NB_WDIR",
   "language": "",
   "name": "rik_ssh_shkpwagn_hdata2_hlrn_de_python_py3_std_gfs2workshkpwagnnb_wdir"
  },
  "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}