Commit 2ca16923 authored by Patrick Wagner's avatar Patrick Wagner

Initial commit

parents
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"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": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"## Load data\n",
"figdir=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/FIGURES/\"\n",
"outfile=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/DATA/GB_arrival.nc\"\n",
"\n",
"ts=[]\n",
"ts_clim=[]\n",
"\n",
"year1=1960\n",
"year2=1961#2007\n",
"\n",
"\n",
"\n",
"for year in np.arange(year1,year2+1):\n",
" dir=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/DATA/GS-\"+str(year)+\"/ariane_trajectories_qualitative_GB_initt.nc\"\n",
" data=xr.open_dataset(dir)\n",
" print(data)\n",
" dummy=((year-1948)*365+(data['ts_gb']+data['init_t'])*5)\n",
" ts=np.append(ts,dummy.values)\n",
" dummy=(data['ts_gb']+data['init_t'])*5\n",
" dummy[dummy>730]=dummy[dummy>730]-730\n",
" dummy[dummy>365]=dummy[dummy>365]-365\n",
" ts_clim=np.append(ts_clim,dummy.values)\n",
" print(ts_clim.shape)\n",
" \n",
"\n",
"year=(ts/365).astype(int)+1948\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEDCAYAAAAWUyJmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm4XVV9//H3J8xjEhRyK0NSAQFFRJTBoRK0YoBqqC1K\nKGW0pQ+gVNQK1JoA/ipQUUCr6M8IgYKAKIOKEBmCioY5zEOQyQC5iCRhEgrk2z/WOsnOyT33nnvv\nPvPn9Tz3ufuss/dee+8zfM9aa+/vVkRgZmZWljGt3gAzM+suDixmZlYqBxYzMyuVA4uZmZXKgcXM\nzErlwGJmZqVqaGCRNFNSv6Q7C2XjJc2W9ICkqySNLTx3hqT5kuZJ2r5QfqCkB/MyBzRym83MbHQa\n3WI5C/hIVdkxwNURsRVwLXAsgKQ9gM0jYkvgMODMXD4e+DKwI7AzML0YjMzMrL00NLBExG+ARVXF\nU4FZeXpWflwpPycvdyMwVtIEUmCaHRFLImIxMBuY0sjtNjOzkWvFGMtGEdEPEBELgY1y+cbAHwrz\nLchl1eVP5DIzM2tD7TR4rwEexwDl5HIzM2tDq7agzn5JEyKiX1If8HQuXwBsWphvE+DJXD65qvy6\ngVYsyQHHzGwEImKgH/Ej0owWi1ix1XE5cFCePgi4rFB+AICkXYDFucvsKuDDksbmgfwP57IBRYT/\nIpg+fXrLt6Fd/nwsfCx8LAb/K1tDWyySzie1Nt4g6XFgOnAS8CNJhwCPA/sARMQVkvaU9BDwInBw\nLl8k6UTgFlIX2PGRBvHNzKwNNTSwRMR+NZ766xrzH1mj/Gzg7HK2yszMGqmdBu+tRJMnT271JrSN\nXj8WfX2TkIQkjj/++GXTkujrm9TqzWuZXn9fNJIa0b/WKpKim/bHrAxS5QTLAZ9tSB+7dRZJRIcN\n3puZWQ9xYDEzs1I5sJiZWakcWMzMrFQOLGZmVioHFjMzK5UDi5mZlcqBxczMSuXAYmZmpWpZYJH0\nWUl3S7pT0nmSVpc0SdLcfG/7H0paNc+7uqQLJM2X9DtJm7Vqu83MbHAtCSyS3gR8GtghIrYjJcOc\nBpwMnBoRWwGLgUPzIocCz0bElsBpwCnN32ozM6tHK7vCVgHWya2StUg39doN+HF+fhawd56emh8D\nXAx8qInbaWZmw9CSwBIRTwKnku7H8gSwBLiNdHOvpXm2yj3voXDf+4h4HVgsaYOmbrSZmdWlVV1h\n40itkInAm4B1gD0GmLWSdrU66+Zg6VrNzKyFWnHPe0g3+no4Ip4FkHQJ8F5gnKQxudVSuec9pNbL\npsCTklYB1o+IRQOteMaMGcumJ0+e7HsumJlVmTNnDnPmzGnY+ltyPxZJOwEzgR2BV4CzgJuBDwA/\niYgLJX0HuCMizpR0OLBtRBwuaV9g74jYd4D1+n4sZlV8PxYbStn3Y2nZjb4kTQf2BV4Fbgc+RWql\nXACMz2X7R8SrktYAzgXeCfwJ2DciHh1gnQ4sZlUcWGwoXRNYGsGBxWxlDiw2FN9B0szM2poDi5mZ\nlcqBxczMSuXAYmZmpXJgMTOzUjmwmJlZqRxYzMysVA4sZmZWKgcWMzMrlQOLNU1f3yQkDfjX1zep\n1ZtnBX6tbDSc0sWaxqlFWmMkx92vVW9xShczM2trLQssksZK+pGk+yTdI2lnSeMlzZb0gKSrJI0t\nzH+GpPmS5knavlXbbdaOBuu6Mmu2VrZYTgeuiIhtgHcA9wPHAFdHxFbAtcCxAJL2ADaPiC2Bw4Az\nW7PJZu2pv/8xUtfVQH/NUyvAeVymt9Q1xiJpHeDPEbFU0luArYFfRMSrI6pUWg+YFxGbV5XfD+wa\nEf2S+oDrImIbSWfm6QvzfPcBkyOiv2p5j7G0MffbN85Qx7ZZYyy1l/Pr285aNcbyK2BNSRsD1wAH\nA2ePot43A89IOkvSbZK+J2ltYEIlWETEQmCjPP/GwB8Kyz+Ry8xsVNZw95mVrt573isiXpJ0KPDN\niDhF0u2jrHcH4IiIuEXSN0jdYIP95Ko24Ly+573ZcLxCrRZGudaoGbAmTJjIwoWPllyfDaYt7nmf\ng8jhwDeAQyPiHkl3RcTbR1SpNAH4XUS8OT9+PymwbE7u4hqiK2xZl1nVet0V1sbcFdY4o+kKqx1Y\nyu0K82vfvlrVFXYUaSD9khxU3gxcN9JKc0D4Qx6vAfgQcA9wOXBQLjsIuCxPXw4cACBpF2BxdVAx\nM7P2UG9X2ISI+FjlQUQ8LOnXo6z7M8B5klYDHiaN26wCXCTpEOBxYJ9c3xWS9pT0EPBintfMzNpQ\nvV1ht0XEDkOVtZq7wtpbL3aF9fVNyqcCr6jscQV3hdlolN0VNmiLJV8/siewsaQzCk+tD7xW1kaY\ndavl15dUl/vMK+teQ3WFPQncAnwMuLVQ/jzw2UZtlJmZda56u8JWG+nFkM3krrD21otdYc26YNBd\nYTYaTe0KK5gk6avAW4E1K4WV04XNzMwq6j3d+CzgO6Rxld2Ac4BzG7VRZmbWueoNLGtFxDWkrrPH\nImIG8MHGbZaZmXWqervCXpY0Bpgv6UhSrq6NhljGzEag1inKAGPGrM3SpS81eYvMhqfewfsdgfuA\nccCJwFjglIiY29jNGx4P3rc3D96v8MwIBsDTciMZHPfgvQ2m7MH7Yd2aOKe7j4h4oawNKJMDS3tz\nYFnhGQeWOtZnzdGSXGGS3p4TUd4N3CPpVknblrURZmbWPeodvP8ucHRETIyIicDngO81brPMzKxT\n1RtY1omIZdmMI2IOsM5oK5c0Jt/o6/L8eJKkufme9z+UtGouX13SBfme97+TtNlo6h3s/uC+haoN\nl+83b7aiegPLw5L+I3/xT5L0JeCREuo/Cri38Phk4NR8z/vFwKG5/FDg2XzP+9OAU0ZT6WD3B691\nNo5ZLe1yv3mzdlFvYDkE2BD4CXBJnh5V6npJm5ASXH6/UPxB4Md5ehawd56emh8DXEy6f4uZtcTA\ntzN2C80q6rqOJSIWke6fUqZvAF8gnbqMpDcAiyJiaX5+Acvva7/snvcR8bqkxZI2iIhnS94mMxtS\nrdsZQ/m3NLZONFTa/J8ySHu+ePOv4ZC0F9AfEfMkTa4Us/K7MgrPrbCKWtvle96bWaca7OLYwe7h\nM9z7/rT0nveSdh1s4Yi4fkSVSv8J7E/KPbYWsB5wKbA70BcRS/MtiKdHxB6SrszTN0paBXgqIla6\n8r/e61h68XqKdtCtx32k15108nUszboupteM9DMy2izazc5u/EhEPF5WZRURcRxwHCwLXp+LiP0l\nXUi6HfGFwIGseM/7A4Eb8/PXlr1NZs21hsck2kyz7vbZC4YavL+0MiHpx4PNWJJjgKMlPQhsAMzM\n5TOBN0qaD/xrns+sg1XGKXwmWbuodXZff//CUi9P6IXT04fqCrs9It5ZPd2u3BVWrpH299bSrce9\nF7uuurErrFnpaEZzU7ZO6QobqsUSNaZ7Uq9dWOnrfazVeu0z1y2GarG8DrxICqNrAZV83SIlo1y/\n4Vs4DI1usXTrL+5ayt7fstdXdotqpHUl7d0iaO/1rUnqGqylOZ+57myx1D62xc9IS7MbtzsHlnKN\ndH8H/xJu30A1mrra+4u7G9eXnnNgGWqd9a2vJdmNzYajdhdae+uFQVWzZnBgaVO1vuTcr9w4zvll\nw1c7vU0vf1bdFbbys23RFTbaPtPGbgM0slne6O0b6frau2uo19aXnmuHrrCRjR11d1dYvfe8t7ZR\nO09Tf38zu2x8gZ/Z4Gp9Vkf6uemcz5y7wmyEOvcCv1rdjGbtrXM+cw4sK+nklODu761Hp55cYKPn\nEzSaw4FlJZ3zq2Bltbe9ky9odAvDyuITNJrDYyzW9pZ/GVRzcDFrRy1psUjaRNK1ku6VdJekz+Ty\n8ZJmK93z/ipJYwvLnKF0z/t5krZvxXZ3o+7tGujkLk2zztaqrrDXgKMj4q3Ae4AjJG1Nylp8db7n\n/bXAsQCS9gA2z/e8Pww4szWb3X26t2ugk7s0zTpbSwJLRCyMiHl5+gXgPmATVry3/az8mPz/nDz/\njcBYSROautFDGvgX8mCD5t3bWjCzXtbyMRZJk4DtgbnAhIjohxR8JFXuErnsnvfZE7msv3lbOpSB\nz1kf7NqS2mMH4PEDM+tULQ0sktYFLgaOiogXJA3nW9Z9GsPSORdXmVlna1lgkbQqKaicGxGVWxD3\nS5oQEf2S+oCnc/kCYNPC4psATw603hkzZiybnjx5MpMnTy55yztV2VcBl82Bz4Zr4PfMmDFrs3Tp\nSwPMb8vNAVb8vixTy3KFSToHeCYiji6UnQw8GxEnSzoGGBcRx0jaEzgiIvaStAtwWkTsMsA6S8kV\nVnaOpJHmHuvkHE6due0+Fp2zvmbW1e7rG3ldjcoV1pLAIul9wK+Au1h+qs5xwE3ARaTWyePAPhGx\nOC/zLWAK6cZjB0fEbQOstw0Dy8hvYtRrb/LuW18z6+q19TWzrnZf38jr6qrA0ijtGVh6bX3NrKvd\n19fMunptfc2sq93XN/K6fKMvMzPrCA4sZmZWKgcWMzMrlQOLmZmVquVX3jfKDTfcwFVXzW71ZpiZ\n9ZyuDSzHHPNVfvOb9YGtqp5Z2orNMTPrGV0bWJL9gL+pKlsKnNCCbTEz6w0eYzEzs1I5sJiZWakc\nWMzMrFQOLGZmVioHFjMzK1VHBRZJUyTdL+lBSV9s9fa0tzmt3oA2MqfVG9BG5rR6A9rInFZvQNfq\nmMAiaQzwLeAjwNuAaZK2bu1WtbM5rd6ANjKn1RvQRua0egPayJxWb0DX6pjAAuwEzI+IxyLiVeAC\nYGqLt8nMzKp00gWSGwN/KDxeQAo2g7gEeLCqzFfem5k1Usfc6EvS3wO7R8Q/58f7AztGxFGFeTpj\nZ8zM2kyZN/rqpBbLAmCzwuNNgCeLM5R5YMzMbGQ6aYzlZmALSRMlrQ7sC1ze4m0yM7MqHdNiiYjX\nJR0JzCYFxJkRcV+LN8vMzKp0zBiLmZl1hrbvCpM0U1K/pDsLZdtJ+q2kOyRdJmndAZ67Oz+/ei7f\nQdKd+eLK01qxL6M1nGMhaT9Jt0u6Lf9/XdJ2+bl39dixWFXS2Xmf75F0TGGZjr/odpjHYjVJP8jH\n4nZJuxaW6ejPiKRNJF0r6V5Jd0n6TC4fL2m2pAckXSVpbGGZMyTNlzRP0vaF8gPzcXhA0gGt2J/R\nGO6xkLRVfr+8LOnoqnUN/zMSEW39B7wf2B64s1B2E/D+PH0QcEKeXgW4A9g2Px7P8lbZjcBOefoK\n4COt3rdGHouq5bYFHio87qljAUwDzs/TawGPkE4EGQM8BEwEVgPmAVu3et8afCwOJ3UjA2wI3NIt\n7wugD9g+T68LPABsDZwM/Fsu/yJwUp7eA/h5nt4ZmJunxwO/B8YC4yrTrd6/Bh+LDYF3AScCRxfW\nM6LPSNu3WCLiN8CiquK35HKAq4G/y9O7A3dExN152UUREZL6gPUi4qY83znA3g3e9NIN81gUTQN+\nCNCjxyKAdSStAqwNvAI8R5dcdFvnsfh4nn4rcE1e7o/AYknv7ob3RUQsjIh5efoF4D7S2aNTgVl5\ntlksf42nkvaTiLgRGCtpAim7x+yIWBIRi0njulOatiMlGMax2DvP88eIuBV4rWpVI/qMtH1gqeFu\nSR/N058gHTCAtwBIulLSLZK+kMs3Jp2uXLEgl3WDWsei6JPkwEJvHouLgZeAp4BHga/lL4yBLrrt\n1mOxaZ6+A5gqaRVJf0n6lbopXfa+kDSJ1IqbC0yIiH5IX7jARnm2Wq9/dfkTdO+x2HCIxUf0GenU\nwHIIcKSkm4F1gP/N5asC7yP9Qv8r4G8l7QYMdH1Lt5y1UOtYACBpJ+DFiLi3UjTAOrr9WOxM+iXW\nB7wZ+Hz+sPXisfgB6YvyZuDrwA2kY9M1xyKPJ10MHJV/rdfaj+p9Vp63F49FzVUMUDbkOjrmdOOi\niHiQ1FxF0pbAXvmpBcD1EbEoP3cFsANwHst/scEAF1d2qkGORcW+LG+tQDpGvXYspgFXRsRS4I+S\nbgDeTR0X3XaqWsciIl4Hlg3O5mMxH1hMF7wvJK1K+iI9NyIuy8X9kiZERH/u8ns6l9f6LCwAJleV\nX9fQDW+AYR6LWkb0GemUFosoRE5JG+b/Y4AvAWfmp64CtpO0Zj6ouwL35Cbfc5J2kiTgAOAyOlO9\nx4K8r/uQ+kWBZc3fXjkW38lPPQ58MD+3DrALqc+5my66ret9IWktSWvn6Q8Dr0bE/V30vvgBcG9E\nnF4ou5x0AgP5/2WF8gMAJO0CLM7dRFcBH5Y0VtJ44MO5rNMMdSwOZODXuNhKGdlnpNVnL9RxdsP5\npAj5CukL4mDgM6SzHO4H/rNq/v2Au4E7ga8Wyt8F3EX6dXZ6q/erScdiV+C3A6ynp44FqSvoovy+\nuJsVz3qZkpeZDxzT6v1qwrGYmMvuIQ1Kb9ot7wtSN/jrpDOXbgduy6/vBqQTGB4AfgmMKyzzLdJZ\nT3cAOxTKD8rH4UHggFbvW6OPBTCBNJayGHg2v4/Wzc8N+zPiCyTNzKxUndIVZmZmHcKBxczMStXQ\nwCJpDUk35tQRd0mansvPkvSwlqcc2a6wTFemWDAz6xUNPd04Il6RtFtEvJSver5B0pX56c9HxE+K\n80vaA9g8IraUtDPpTJZd8pkZXyadOizgVkmXRcSSRm6/mZkNX8O7wiLipTy5BimQVe4NPNCFN12b\nYsHMrFc0PLBIGiPpdmAh8MuIuDk/9ZXc3XWqpNVyWU+kWDAz62YNv/I+0tXO75S0PnCJpLeSzoXu\nzwHl/5OybH6FUaZYkO95b2Y2IlHird2bdlZYRDwHXA9MieVJ0F4FziJl0ITBUyzUlVZguBcSTZ8+\nveEXK7mO9qunm/bFx6w962hmPaP9K1ujzwp7Y+FGMmsBfw3cn3PUVFKO7E26Ghq6P8WCmVnXa3RX\n2F8As3K+ojHAhRFxhaRrJL2R1MU1D/gXgPzcnpIeAl4kpaYgIhZJOhG4hdQFdnykQXwzM2szjT7d\n+C7SKcLV5R8aZJkja5SfDZxd1rZVTJ48uexVuo4OqKeb9qVZ9Xhf2reedtNVucIkRTftTzfo65tE\nf/9jDa9nwoSJLFz4aMPrMetGkogSB+8dWKyh0jBaM14TNWQQ0qwXlB1YnCvMusQaSGroX1/fpFbv\npFlHcIvFGqqZLZbG1+NWkXUnt1jMzKyttSq78SRJc3Om4h/m2wgjaXVJF+Tsxr+TtFlhXcfm8vsk\n7d7I7TYzs5FraGCJiFeA3SLincD2wB45a/HJwKkRsRXpVpiH5kUOBZ6NiC2B04BTAHIamE8A2wB7\nAN/OF1eamVmbaUV24wB2A36cy2eRrr6HlN14Vp6+GPhgnv4YcEFEvBYRj5LuvVxJA2NmZm2k6dmN\ngd+TUrVU0udXMhhDIYtxRLwOLJG0Ac5ubGbWMZrRYlmau8I2IbUythlotvy/VhbjurIbm5lZ6zU8\nbX5FRDwn6XpgF2CcpDG51VLMVFzJbvxkvuPk2JwnrFbW45XMmDFj2fTkyZN7NqWCmVktc+bMYc6c\nOQ1bf0OvY8mJJl+NiCU5u/FVwEnAgcBPIuJCSd8B7oiIMyUdDmwbEYdL2hfYOyL2zYP35wE7k7rA\nfglsWX3Riq9jaT++jsWs/ZV9HUurshvfB1yQMxbfDszM888EzpU0H/gTsC9ARNwr6SLgXuBV4HBH\nEDOz9uQr762h3GIxa3++8t7MzNqaA4uZmZXKgcXMzErlwGJmZqVyYDEzs1I5sJiZWakanTZ/E0nX\nSro3p83/dC6fLmmBpNvy35TCMgOmx5c0RdL9kh6U9MVGbreZmY1co6+87wP6ImKepHWBW0kZjD8J\nPB8RX6+afxvgfGBHUtqWq4EtSRcpPAh8iJTK5WZg34i4v2p5X8fSZnwdi1n766gr7yNiISmrMRHx\nQr7ivpKVeKCdmEpOjw88mq/A3ynPOz8iHgOQdEGe9/4B1mFmZi3UtDEWSZNIN/u6MRcdIWmepO9L\nGpvLaqXHry4vpto3M7M20pTsxrkb7GLgqNxy+TZwQkSEpK8ApwKfonZ6/IEC4IB9Es5ubGY2uI7O\nbgyQ72f/M+AXEXH6AM9PBH4aEdtJOgaIiDg5P3clMJ0UcGZExJRcvsJ8hXV5jKXNeIzFrP11Yq6w\nHwD3FoNKHtSv+Dhwd56+HNhX0uqS/hLYAriJNFi/haSJklYnZT2+vAnbbmZmw9TQrjBJ7wP+Abgr\n3544gOOA/SRtDywFHgUOg0HT478u6UhgNikYzoyI+xq57WZmNjJOm28N5a4ws/bXiV1hZmbWQxxY\nelhf3yQkNfTPzHqPu8J6WHO6qdwVZtbu3BVmZmZtzYHFzMxKVVdgkbSOpDF5+i2SPiZptTqWq85u\n/JlcPl7SbEkPSLqqkNIFSWfk7Mbz8inJlfIDc2bjByQdMPxdNTOzZqhrjEXSrcBfAeOBucAtwEsR\n8Q9DLFcru/HBwJ8i4pScAn98RBwjaQ/gyIjYS9LOwOkRsYuk8bnOHUid6bcCO0TEkqr6PMYyDB5j\nGX4dfn9ZN2rVGIsi4iXSVfLfjIi/Bd461EIRsTAi5uXpF4D7SOnwpwKz8myz8mPy/3Py/DcCYyVN\nAD4CzI6IJRGxmHSh5LJ7uJiZWfuoO7BIeg/pKvqf57JhXbVfyG48F5gQEf2wLLX+Rnm2WlmMa2U9\nNjOzNlNvYDkKOBa4JCLukfRm4Lp6K6nObkztPovqplilf6NW1mMzM2sz9bY6JkTExyoPIuJhSb+u\nZ8Gc3fhi4NyIuCwX90uaEBH9eRzm6Vy+ANi0sPgmpDtGLgAmV5UPGNicNt/MbHBtkTZf0m0RscNQ\nZTWWPQd4JiKOLpSdDDwbESfnFPjj8uD9nsARefB+F+C0AQbvx+Tpd+XxlmJdHrwfBg/eD78Ov7+s\nGzX11sT5LK09gY0lnVF4an3gtaFWPkh245OBiyQdAjwO7AMQEVdI2lPSQ8CLpLPHiIhFkk4kBZQA\njq8OKmZm1h4GbbFIegdpwP0E4MuFp54HrouIRY3dvOHplhZLX98k+vsfa1JtbrEMp45ueH+ZVSu7\nxVJvV9hqEfFqWZU2SrcElm5LNd9N+9IN7y+zak3tCiuYJOmrpGtX1qwURsSby9oQMzPrDvWebnwW\n8B3SuMpupIsYz23URpmZWeeqN7CsFRHXkLrOHouIGcAHG7dZZmbWqertCns5J6Gcn+89/wTLr5Y3\nMzNbpt7B+x1Jeb7GAScCY4FTImJuYzdveDx4P+yamlBPd+1LN7y/zKq1JAllRNycU7EsAT4dER+v\nJ6hImimpX9KdhbLpkhZIui3/TSk8d2xOmX+fpN0L5VMk3Z/T5n9xeLtoZmbNVO/9WN6eL3C8G7hH\n0q2Stq1j0bNImYmrfT0idsh/V+Y6tgE+AWwD7AF8W8kY4Ft5PW8Dpknaup7tNjOz5qt38P67wNER\nMTEiJgKfA7431EIR8RtgoIsoB2pyTQUuiIjXIuJRYD6wU/6bn08aeBW4gOVp9s3MrM3UG1jWiYhl\nSR8jYg6wzijqPSLfIfL7hbtH1kqNXyuVvpmZtaF6zwp7WNJ/sPzalf2BR0ZY57eBEyIiJH0FOBX4\nFLVT4w8U/GqOoDq7sZnZ4Nolu/F44Hjg/aQA8CtgRj25wiRNBH4aEdsN9lzOchwRcXJ+7kpgeq5v\nRkRMyeUrzFe1Pp8VNryamlBPd+1LN7y/zKq1JKVLDiCfGWEdotAakdSX7xoJ6VbHd+fpy4HzJH2D\n1NW1BXATqcWyRQ5CTwH7AtNqVfZ3f3fgCDezfu9+9/Yce+xnG16PmVknGiq78U8Z5Gdg8eZfNZY/\nn3SDrjcA/aQWyG6kjMlLgUeBwyq3KZZ0LHAo8CrpbpOzc/kU4HRSkJkZESfVqC/g7ME2qQSLGDv2\nFBYvfrJhNbjF0q71uMXSy5qVdXzChIksXPhow+spamp2Y0m7DrZwRFxf1oaUIQWWRn/wn0TakoiX\nGlxP93wZd9O+OLD0rmb+4Gv2+6zZXWGPRMTjZVXWLVJQaeQLX9rra2bWdEOdbnxpZULSjxu8LWZm\n1gWGCizFn86+94qZmQ1pqMASNabNzMwGNNQYyzskPUdquayVp8mPIyLWb+jWmZlZxxm0xRIRq0TE\n+hGxXkSsmqcrj4cMKjWyG4+XNFvSA5KuKqR0QdIZObvxPEnbF8oPzJmNH5B0wEh31szMGq/eXGEj\nNVB242OAqyNiK+Ba4FgASXsAm0fElsBhwJm5fDzwZWBHYGdgejEYmZlZe2loYKmR3XgqMCtPz2J5\npuKpwDl5uRuBsZImkALT7IhYEhGLgdnAFMzMrC01usUykI0qV9rn1C6VWxzXymJcK+uxmZm1oVYE\nllqqrwqsXOZaK+uxmZm1oXrT5pepX9KEiOiX1Ac8ncsXAJsW5tsEeDKXT64qv46aZhSmJ1ctamZm\nbZE2f1QVSJNIqfHfnh+fDDwbESfnFPjjIuIYSXsCR0TEXpJ2AU6LiF3y4P0twA6kFtYtwLvyeEt1\nXU3JFZZ64hqd0qV78mt10744V1jvcq6w+jW0xVLMbizpcVJ245OAH0k6BHgc2AcgIq6QtKekh4AX\ngYNz+SJJJ5ICSgDHDxRUzKx3NSvzsNWn4S2WZnKLpR3r6a596abPSzdpTmvCLZZ6tdPgvZmZdQEH\nFjMzK5UDi5mZlcqBxaxuayCp4X99fZNavaOl6eub1JRjZu3Fg/fD5sH79qujWfV07+BtozTzFF2/\n/qOo0YP3ZmbWzloWWCQ9KukOSbdLuimXDTulvpmZtZdWtliWApMj4p0RsVMuG1ZKfTMzaz+tDCwa\noP7hptQ3sxFqxsC69aZWBpYArpJ0s6RP5bIJdabUd+p8s1FKKVCiwX/Wi1qR3bjivRGxUNKGwGxJ\nD1D7nTiM1PkzCtOTcXZjM7MVdXx247o2QpoOvAB8ijTuUkmpf11EbCPpzDx9YZ7/fmDXSuumsB6f\nbtx29XhfRlJPMz6X3ZZfq5v2xacbj4CktSWtm6fXAXYH7gIuBw7Ksx0EXJanLwcOyPPvAiyuDipm\nZtYeWtV8KUTLAAALiElEQVQVNgG4JLUwWBU4LyJmS7oFuKjelPpmZtZ+2qIrrCzuCmvHerwvw7cm\n8EoT6oHuOWbd9Pp3fldYKwfvzWxAr9C8L2Oz8jmli5mZlcqBxczMSuXAYmZmpXJgMTOzUnVUYJE0\nRdL9kh6U9MVWb4+Zma2sYwKLpDHAt4CPAG8DpknaevRrnjP6VbiODqynGXV0Wz3NqKNZ9TSjjmbW\n0146JrAAOwHzI+KxiHgVuIDl2Y9HYc7oV+E6OrCeZtTRbfU0o45m1dOMOppZT3vppMBSneF4Ac5w\nbGbWdjrpAsm6Mhyvv/5Hh7XSl19+gDXXvLXu+SNe5vnnh1WFmVlP6ZiULjn55IyImJIfHwNERJxc\nmKczdsbMrM2UmdKlkwLLKsADwIeAp4CbgGkRcV9LN8zMzFbQMV1hEfG6pCOB2aSxoZkOKmZm7adj\nWixmZtYZOumssK4g6bp8vxnrQZLOknRCq7fDrJG6PrA08otc0qOSXpL0nKTn8/++BtTxsqQNqsrn\nSVoqabOS65sj6VlJq5W83mbvR1MDeIPfZw15TarqeL+kGyQtlvSMpF9LeleD6jpI0p2SXpT0pKRv\nSxpbx3JLJb15iHkelbRQ0lqFskMlXVfGtlfV85KkJfm1+Y2kw5Tu99zzuj6wNFgAe0XE+hGxXv6/\nsAF1PAJMqxRI2pZ0N6hh92PmkyBqPTcReD+wFPjYIPON5H1T6n70inpfk1HWsR7wU+B0YDzp+rDj\nacDdxiR9Dvgq8DlgfWAXYCLwS0lDjfnW8z4JYBXgX0ew7HBUPvtjSdt/EvBFYGbJ9XSkngksksZJ\n+qmkpyX9KU9vXHj+Okkn5F8ez0m6svrXda1VD1DXLvnX3yJJt0vatWqWLSTdmH8dXiJp3BB1nAsc\nWHh8IDCrUN+ekm7Lv54ekzS98NzE/EvvEEmPAdcMUs8BwO+As4GDCus4K/+q/Lmk54HJuew7kmbn\n43VdHa2O0ezHzyQdUVyZpDskDfplK+lASb+uKlv2yzfvx7fy+p+T9DtJfznEfoy4vhGo9Zqs0EKq\nrlfS7kp59RZJ+u/c6qnVonoL6dT9iyJ5JSKujoi787oOkXRv/tz8ovg65337tKTf58/WKbV2JAew\nGcCREfHLiHg9Ih4HPkH6ct5f0hhJx0l6KL8eN0vaRNL1pM/anbl8n0GO2X8Bn5O0/gDb8F5JN+Xj\ncqOk9+TyT0q6uWrez0q6dJB6RDpwz0fEz4BPAgdKequk1SV9Lb+Pn8qfnzUK656avxuWSJovafdB\n6uk4PRNYSPv6A2BTYDPgJVLusaJppC+7DYE1gM8PtxJJbwJ+BpwQEePzOn4s6Q2F2f6R9CXxF8Dr\nwDeHWO1cYD1JW+XWwieA/2F5UHsB+Mf862kv4F8G+ML9ALA1KddaLQfk9Z4PfETShoXnpgEnRsR6\nwA25bD/SL9s3AHcA5zVwP2aRjhsAkt4BvAm4Yog6YeVfq9WP9wWmA+OA3wP/r451jqa+4RjsNRmw\nXklvBH5E+gX9BtJp+u8ZZLkHgdclna2U6HXZDx1JewPHAHuTPhe/Bn5YtfzewA75b+ogAey9pM/V\nJStsdMSLwC+ADwNHk76gp0TE+sAhwIsRUflx9vbcM/CjQfbnFlIulS8UCyWNJ302TyMdl28AP8/l\nlwNvkbR5YZFpDP2eLu7HzaSMIH8FnAxsAWyX/28MfDlvx06k9/Pn8nv9A8Cj9dbTCXomsETEsxFx\nSf419iKpOf6BqtnOiojfR8QrwEXA9nWs+lKlPtZnJf0E2B/4eURcleu9hvRG37OwzLkRcV9E/Bn4\nD2Afaci+2cqv/Q8D9wNPFvbtVxFxT56+m5RHrdhKCmB6RPw579tKJL2fFHAviojbgIdIgaPisoiY\nm+uorOPnEXFDzt3278B7VGgFlrwfl5FaepUP/v7AhRHx2hD1Dbi7VY9/EhG3RsRS0hdJPa/7aOqr\nb6GhX5Na9gDujojLImJpRJwB9NeaOSKeZ3l32/eAP0q6VNJGwD8DX42IB/PxOQnYXtKmhVWcFBFL\nImIB6Ut7GgN7I/BMXk+1p0iB61PAlyLiobxtd0XEosJ89R7L6cCRVT/o9gIejIjz83G5gPQe/Gj+\nLF5e2XZJWwJbkboIh+NJUtD6J+Cz+bi8SDpuleNyCOlyiWvzPj4VEQ8Os5621jOBRdJakr6rNOi2\nGLgeGFf1hV4cH3kJWLeOVU+NiA3y38dJTfpPFILNIuB9QHFQv5jz7DFgddKHbjD/Q/pSOQg4p2rf\ndpZ0be6KWAwcNsD6Fgyx/gOA2YUP8Q9ZsdvqDysvsrwsf3ieJbUiSt+PiPhfUrDfP79m00hBqgwj\ned2bYajXpJY3sfLrNejrHxEPRMQhEbEZKXv4m0hBYiJweuX9DPyJ9EOl+AOiuO7HqP0eeAZ4owYe\no/uL/PwmpFbjqOQfKD8Djs1Fytv1WNWsj7F8X85n+Zf/fsClEfHyMKvemDTGszZwa+G4/YIUcCD1\nmox6H9tZzwQW0mDhlsCOETGO5a2V0Z7FUb38H4BzCsFmfB7Y/6/CPMVfexOB/yV9qGrKfdGPkH6N\n/qRSnP+fB1wKbJz37bsDbFfN7hhJa5K6pXbN/cFPAZ8F3iFpu0GW37SwjnWBDSi0QBqwH+eQWiof\nInWP3DhYXdmLpA95ZTtLPWuvUfXV8ZqsUA8r/nB5ihXfY5C+sOuSfz3PArYFHgcOq3o/r1tpvWbF\nujaj9nvgd6QTAj5eLJS0Dun9cDXp87P5youOyAxSy2Fj0nvsCWBS1Tyb5XJIF1+/MXez7ksKNHWT\ntCMpeF1K+oHytsJxG5e7vaDcfWxLvRRY1gP+DDynNCg/o0H1/A/w0Tx4OkbSmpJ2zWMvFftL2lrS\n2qQxih9FfVeqHgJ8MDfbYfmX7rrAooh4NfffVneXDBU8/xZ4DdgGeEf+25rUn37AIMvtmQdDVwdO\nBOZGxBODzD+q/chfZkuBU6m/tXIH8DZJ2+XB0+k09iy0suob6jW5Hfh4bolvARxaWPbnwLaSPiZp\nFaWMFRNqVZTHvI6udGPmbq5ppEBwJnCcpLfm58ZK+vuqVXxB6eSYTYGjSF2YK4mI54ATgG9K+oik\nVSVNIrVEHye9pjOBE/M+IenteQwEUsuy7pMgIuL3wIXAZ3LRL4AtJe2bj8snScf3Z3n+14GLSYP/\n44Ff1lOPpPUk/Q2pRXluRNwFfB84rTImJmnjwgD9TOBgSbspeZOkrerdr07QK4ElSAN1a5NaBr9l\n5UHfkXz4V1om9zNPBY4D/khqan+e5cc6SB+gWaRfdquTPoxD1hERj+S+9urnDid9GJcAXyJ9mAbd\nzioHAD+IiCci4unKH/DfpC/3Wqcon08K0H8C3gn8Q4P3A1KrZVtSAB9KRMR80pfZNaRB6l8Pvsio\nlFnfUK/JN4BXSV+2Z1E4HhHxJ2Af0hfkM6SAdAu1Tx9+HtgZuFHprL/fAneSBpcvI40PXJC7J+8E\nplQtfxlwK3AbaUziB7V2KrfcjwO+BiwhBa/HgL/OY3VfJwWa2fl98H2gck3K8cA5uXupOrgtq6Lq\n8Qmkz31ExLPA35A+j8/k/3vl8oofklrEF9UYCyr6ad7Gx0ldbl8j/WgC+DfSmNjcfNxmk86+qwzy\nH0zqalxCOtGg1Ou4Wq3rU7pIuhU4PiIub/W2dBNJZwF/iIgvN7nefwT+KSKqT7yonq+pr3s7v8/y\nmNQCYL+IuL7kdS8FtoiIh8tcr3W2rm6xSHob6dfa7a3eFhu93HV4OGnsZbD5mvq6t+P7LHfFjs3d\ncf+ei+cOtoxZWbo2sEg6CbgS+LeIGOiMJhudpjZ1c//006SB6errKIrzNfV1b+P32XtIZx49TTrN\ndmqtU81Hqbu7PGxEur4rzMzMmqtrWyxmZtYaDixmZlYqBxYzMyuVA4uZmZXKgcXMzErlwGJmZqX6\nPztliAaA9MxoAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x2aaacc18b278>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## Plot\n",
"\n",
"fig,ax=plt.subplots(2,1)\n",
"out=\"gb_arrival.pdf\"\n",
"\n",
"yebin=np.arange(year.min()-0.5,year.max()+1.5)\n",
"mdays=[0,31,28,31,30,31,30,31,31,30,31,30,31]\n",
"mlabel=['Jan','Feb','Mar','Arp','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec',]\n",
"mobin=np.cumsum(mdays)\n",
"\n",
"nyears, bins, patches=ax[0].hist(year,yebin)\n",
"nmonth, bins, patches=ax[1].hist(ts_clim,mobin)\n",
"\n",
"cbin=bins[1:]-np.asarray(mdays[1:])/2\n",
"ax[1].set_xticks(cbin)\n",
"ax[1].set_xticklabels(mlabel, fontsize=12)\n",
"\n",
"ax[0].set_ylabel('Floats')\n",
"ax[1].set_ylabel('Floats')\n",
"\n",
"plt.savefig(figdir+out)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"## Save to file\n",
"\n",
"dataout=xr.Dataset({'count': (['year'], nyears), 'count_clim': (['month'],nmonth)},\n",
" coords={'year': np.arange(year.min(),year.max()+1)})\n",
"dataout.to_netcdf(path=outfile)"
]
}
],
"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.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"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": 4,
"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": 41,
"metadata": {},
"outputs": [],
"source": [
"## Set directories\n",
"figdir=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/FIGURES/\"\n",
"outfile=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/DATA/GB_transit_days_median.nc\"\n",
"dir=\"/gfs1/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/DATA/\"\n",
"\n",
"\n",
"## Set parameter\n",
"year1=1960\n",
"year2=2007\n",
"\n",
"\n",
"## Calculate median\n",
"i=0\n",
"median=np.zeros(year2-year1+1)\n",
"uquartile=np.zeros(year2-year1+1)\n",
"lquartile=np.zeros(year2-year1+1)\n",
"for year in range(year1,year2+1):\n",
" file=dir+\"GS-\"+str(year)+\"/ariane_trajectories_qualitative_GB_initt.nc\"\n",
" data=xr.open_dataset(file)\n",
" tsgb=data['ts_gb']\n",
" median[i]=tsgb.median()\n",
" lquartile[i]=np.percentile(tsgb.values,25)\n",
" uquartile[i]=np.percentile(tsgb.values,75)\n",
" i=i+1\n",
"\n",
"median=median*5\n",
"uquartile=uquartile*5\n",
"lquartile=lquartile*5"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD8CAYAAABthzNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsnXd4U9f9uN8j7703NtjYGLP3CiYJZEIW2atNs5NmNF3pbtKV9tv215Hd7JSMZhJIQiYhCcNAGDbEZtjYeMp7b1s6vz+O5CnJsiV5wH2fR4/h3nPvObKl+zmfLaSUaGhoaGhoDEQ31gvQ0NDQ0BifaAJCQ0NDQ8MimoDQ0NDQ0LCIJiA0NDQ0NCyiCQgNDQ0NDYtoAkJDQ0NDwyKagNDQ0NDQsIgmIDQ0NDQ0LKIJCA0NDQ0Ni7iP9QIcITw8XE6ZMmWsl6GhoaExodi/f3+1lDJiqHETWkBMmTKFffv2jfUyNDQ0NCYUQohCe8ZpJiYNDQ0NDYtoAkJDQ0NDwyKagNDQ0NDQsIgmIDQ0NDQ0LKIJCA0NDQ0Ni2gCQkNDQ0PDIpqA0NDQ0NCwiCYgxors96A6b6xXMbHo7oB9L0B351ivREPjtEATEGNBUwW89T3433XQ1TbWq5k47HsRPvgh5Gwa65VoaJwWaAJiLDiyGZBQfRy2/WmsVzMxMBpg95Pq3/nbxnYtGhqnCZqAGAuyN0JEGiy8GXY9DkV7xnpF458j70N9IfhHQf6XIOVYr0hD45RHExCjTaMeCnfBzPVw3h8gKB7euxs6W8d6ZeObjCcgJBHOfBAaS6E6d6xXpKFxyqMJiNHGbF6aeRl4BcClj0PtCfjij2O9spFTXwRbfqp+uoLivVCyF5Z9H6auUcdcaWZqKIEPfwztja6bQ0NjAqAJiNEmeyNEzoSIVPX/pDNh8e3Kvl64a2zXNhKOfABPr4S9z8CXf3HNHBmPg3cwzL8BQhMhZAqccKGA+OY59fr0166bQ0NjAuAyASGESBVCZPZ5NQohHhBChAohPhNC5Jp+hpjGCyHEo0KIPCHEISHEAletbcxoLIOiDGVe6ss5D0PIZHjv+9DZMhYrGz7dHfDRz+GNG5TpZ8alcOhNaCp37jx1J5X/YdHN4OmnjiWdDSd3gKHLuXOB8m1kbwQ3LzjwMuR97vw5NDQmCC4TEFLKY1LKeVLKecBCoBXYCPwc2CqlTAG2mv4PcCGQYnrdATzlqrWNGebwzJmX9T/u5Q+XPgl1BfD570Z/XcOltgCePw/2PAVL74ZbP4U1D4GxG/Y+69y5dj8Fwg2W3Nl7bOrZ0NkEpfudOxeAPlMJpfP/BBHTYdN90Fbv/Hk0NCYAo2ViWgOckFIWApcCL5uOvwyYn5aXAv+Vit1AsBAiZpTWNzpkb4So2RCeMvjclDPUw3bvf6Bg++ivzV6y34P/rFLC7JpX4cK/gLsXhE2F6etg3/PO04La6uHABph1BQT2+ShMSQeEa8xM2RtB567mvOxJaK6AT37l/Hk0NCYAo9VR7lrgddO/o6SUegAppV4IEWk6HgcU97mmxHRMP0prdC0NJVC8B1b/xvqYNb+F3E/g9WvBL3zweZ07rP71YBOVMynZBx89CK01g89Jo3JET1oMV74AwQn9zy+/B45+AFmvw+LbHF/L/pegqwVW3Nv/uG8oxM5Xjuqzf+H4PGbM5qWks9UcvqGw8oew/e8w4xKYdr7z5tJQG4nPfgvu3kpjs5ev/qo+nxf+n+vWZouiPSphs8vKRmjZ92HpnZbPTTBcLiCEEJ7AJcBQ32Rh4digYHchxB0oExQJCQmDLhi39JiXbDzcPX3VrjzjCTBasK/rs5TJI27h4IezoxiNsPsJ+PxhCIiBySssj1v4PVhxP7h5DD6XsBxiF0DGk7DwFtA5oKB2d8Ke/0DimRA9e/D5qWfDjn9BewN4B418nr6UHVAC8Myf9x4780E49hFsvh/u2Q0+Ic6Z63SnIltVE6g+DkKnBLGlTdFAutpg57+hsxnil8Ksy12+1H50NMO7t6nEzSkrB5/XH1LJr/NuUKbjCc5oaBAXAgeklBWm/1cIIWJM2kMMUGk6XgLE97luElA28GZSymeAZwAWLVo0cbKlsjdC9BxlirFF1Ay47AnL5+oK4akVsOle+O4mEJZk6ghorYWNdyntZfpFKvR2JA9CIdRu/+1b4PhHyuQ0UnLeg6YyuPjfls8nnQ3b/59yVjsyT1+yN4LOA6av7T3m7qVMTc+uVk75y//jnLlOV6SEA/9VWqp3EFz4V/XvI5th0S1DX5+3VQkHvwgVijxlJfhHDn2ds/jst1BfDLd8DAnLBp8v3gvPnwuZr8HSO0ZvXS5iNHwQ19FrXgLYDNxk+vdNwKY+x79rimZaBjSYTVETnvoiKPnGcdNQyGSVXFfwlSpa5wwKM1SYav42uPBvcM0rju2S0y5VyX8ZVoScPUgJux6D8FRIPsfymPgl4OGrsqqdgZTKvzJ19eD3HzsPVv0EDv0Pjn7onPlORzqa4J3b4P371cP1rh2w5A4IS1HC2R6yN4JPKHznPWWi+uCHo5dVf2Kb8rEtv8eycAD1uZy0RGnjRsPorMuFuFRACCF8gXOBd/sc/gtwrhAi13TOHDy/BcgH8oBnge+7cm2jirXopZGw8GZIOgs+/Y2KthkpRqPagb+0Dtw8VSTS0jsc10rc3GHpXVC4E0oP2B7b3am+5ANfJ76A8kPqi2jNTOXupcxgznJUl+6HhmLrQjz9JyrA4P0HlMblLLraLP8OTrXMen2WCm7Iflf50W7cqHb+Qqjf+ckd0Fxp+x5dbcrcl3YxRM+C1b9SPq/Db9u3BkcESXsjbL5PCbPVQ+THLL9HfTePbRn5fOMEl5qYpJStQNiAYzWoqKaBYyVwjyvXM2Zkb4SYeRCa5Pi9hIBLHocnl5tMTZuHb+tvroKNd6gH8czLlRnHO9DxtZlZ8F2VNJfxuHJmD0RKpWFs/T0YOizfwzcc5lxje56ks+HTX6kAgKBJjq05e6MSlKkXWj7v7gnrn4JnzlJBBFe+CEFxjs15YANsvtf6+YQVcP4flc9pItNQqsKifULhex8O9m/NXA9f/1WZmWwFN+R+phzDZiG+/F6VI7PlJ5CYDgHRlq+ryIF3b4fINLjiuZG9h09/rUq83PIpePjYHpt2MQRPVnXW0i4e2XzjhNGKYjp9qStUu9NznJjfEBwPFzyidjTfPDc8W2fBdqXmt9XBRf9UGomzfBlmvANh4U0qh+Gc36n1mmmtVbWnjn8M0y6w7gyPXwYe3rbnmXq2+nliGyz4zsjXazSazEtrwCfY+rjo2XD5s0owP70S1v8Hpp03wjkNSoOLSIN51w0+39kC3zyvfB+zr1J5Jn1/jxOJ3E+hux1ufEf52AYSmabMidnv2RYQ2RvBN8wU5gzo3OCyp9Tf4v0H4LrX+3+W+/o7DF1Q8S2s+mlvFQN7yftcJU2e8QDELx56vM4Nlt0NH/9cRQVOWjS8+cYTUsoJ+1q4cKEc9+z4l5QPBUpZW+Dc+xqNUm64XMo/RktZnTf0eEO3lNv+LOXDwVI+ulBK/WHnrmcgdYVSPhwi5ce/7D1WmCHl/5sh5e/Dpdz9tHoPjmA0Svm3FCnfutmx+xTtUX+jzP/ZN77quJRPrlDXfPJrKbs7hz9nzvvq+m/ftT6mrUHKzx6W8g+RUv4+QsrPHlLHJhpvfEfKv0+3/ff+4hEpHwqSsrHc8vmOFvVZ3/yDwed2Pa5+lwdf6z3W3ijl27eq4y9fImV5tvodbr5/eGtvq5fy/6VJ+fgSKTvb7L+uvVHKR+KlfPOm4c03SgD7pB3PWK0Wk6vJ3qhCP0OmOPe+QsDFj6qom033qF2wNZrKYcNl8OWfYfbVcMeXyobrSoITVPmNA/9Voajb/wEvrlU+ils/VXHijmouQih/TP5Xtt//UJhLa1gzLw0kPAVu+1xF3ex6FF68cPiFCjMeh6AEmG7DBOEdCOc8BPfuU7/LHf+ER+fDwVeHN9dYYjRAwddK27P19555GSBNxSwtkPspdLVa9hEtvUuFWH/0M1XORn8I/nMmfPuOyd/xrtJc5l4LWf+Dlmr71//JL9X357Inh9Zo++IVoLTonE3KijBB0QSEvUipbN3DobYAyg66LrEtKE5lMhdlqId/8d7Br2/fVSp48Tdw6ROw/unRi89ecS90NMJTZ8DW36lkszu/VkluziLpLGitVuYDS0gJNSesOyjN5qXkc4bnh/HwUSa6K1+EyqPwdLr9EU4l+9XfbNndSmAORXA8XPEs3P6F8mNt+j5UHbd/rWOJPkuZM5POsj0uMk2Z26xFM2VvVKGtk88YfE7npj7bhk545Qp47hwlTG76QJmUdG5q3PJ7lanrm+ftW/vxT+HgK7DygZH5gZbepXI89kzc0GhNQNjL3mfgnzPVLqXbimO1L3WF8M6t6gMy41LXrWvudTDtQuXke/7cwa+3b1YO3zu+hPk3Ot/fYIu4hTB5pYpOMT9MnZXUZibpLPXTUvnv1lp4/Tp4bIHK8+hoHjymZK/KtxipEJ91Odz1tdIQ/3e92i0PRcbj4BU4fL9J3EK49lWlNTorzNnVmP8uSWcNPXbmelXReGDBx84WpUGkXWJdoIZNhXN/B5U5ymF91w5VvqYvEdMg5Xz45lnoare9lrY6FY4bOQPO/NnQa7dEUJwKAjFr0RMQzUltD4ZuFZfvGwZ7noai3XDVi9ajko58oHZ5UqqHYshk161NCLj6vyqsVFqIuxZuKmZ7qMgLV3HtKypk09GIH2sExqqieie2wRk/6D1etEcl7DVXqLpKh95QwQJXvdTfvNZjXrpg5GsITYKbP1Ka2qZ74O5dysRgifoiZXZY/n3rY2zhH6k2HJmvwZrf9Fa4Ha+c2AZRs+xLZpt5GXz5CORs7h94Ycu81JcldygNI3KG9ci+5ffAfy9Rn4eFN1keA/DxL9TG5rrXVUj1SFl+Dxx+UwmJFfeN/D5jhKZB2MORTSpG/pLH4NrXVKG6/5w5WB3u7lAaxhs3qIfGnV87J/dhKNw9lY03+ZzBr6lnj51wAJV05irhYCbpbGWy6WpXJqMd/1J+AbO/48oX4KbNytz17Gq1+5ZSjc3ZBCnnjuxh3RdPX2Wnri9W2bbWMJsblt418rkW3wYdDcrGPp7pbFX1x5LOsm98RKrqlTLwe5W9EfwirUe8mRFCCX9bYd+Jq1Q0WsYT1s2OR7eoemKrfuK4OTR2noq62v20a8rTuxhNQAyFlCqeOTRJmXKmr1Pqa0SqqiXzwY/Ug6k231QC+2lVlfWWT1RzGw3XM/VsZVs++gG8djV8/hCkXaQEdJyprUjiKrhrpzI7fPBDpV3kfQ5Neuf5iBKWqR3jvhdUjslA2htg/8tqPkfyNhKWqV3y3mfHd2/uol3KL2AOR7aHmeuVsG80VdnpaFa+gBmX9voSHEEI5YuoPma510drLXzwgEqKTP+J4/OB+kw0lvQmzE4gTk8B0d6gCn7ZE/lStFsVcVv2/d6dSXCCMimsuF+l3j9zptIoBpbA1hgdJq9QlW7fuVX5ANb9P7jq5cH+Dv8IuOEdVTU3ZxO8fo2qJDrNAfPSQFb/WmXbbrpvcMvSA/9VfSwGVqcdLkLA4ltVtrkjPTGkVA9JVzWpOrFNJR8mDLHz74s5minHFM2U+wl0tzk30GPm5aog5a7HBp8zVzJe/5TSzJ1ByvkQlqzmG88C3QKnp4A49pEyA+yzI5ohw1S4bt4N/Y+7eai6SNe/qWyVEalKs0i7yDVr1rCOV0Dvl/C2z5UJxpozXqeD9B+rjN6AWJh1pXOjujx8VPJWU1n/lqWGbmVmmLzSOVFcc64BT3+VKDkSWmtVRvgrV6jy2a4g/0tVcdXT1/5rwlPU7t1sZsreCP7R1msfjQR3TxVmXfAVlB/uPZ6zGQ6/pZzSlioIjxSdTvkf9Jnq/hOI01NAzLkGks9VQqI23/q4mhMqdHHRrdY/5NPOhx8fhVs/c34Jbg37uWaDyheImWPf+MnL4YHDcMmjzl9L/GL1QOjbsjTnPWVmWO6kajJeASqu/9t3h18bqmi3CsvN26pyMbI3On9n21ypQo+HY14yM/NSKN4NVcdUeQ1nmZf6svB74OHXW1SypVqZHmPmqtLjzmb+d1QRvy0/hcaJU4P09BQQQqj6QzoPeM9Gktnup5SmsOR22/dz9xrd8FGNwejchv830Omc/+Axc9YvVfkIc8vSjMeVhuNMc9aiW1Utq4Ov2DfeaBycsHjWz6C+UOXrOJP8r9TPpBEIiBkmc9Kme5RvyRV5RD4hKuz78Nvqgb3lJ8r0fNlTlnudOIq5LEh3u/JxTBBT0+kpIKBPktku1eZzIK21kPmqqoNjrQiYhoY1PLyVHbu5Qplxyg7292M5g6gZyr6/7/mh/WnNVfDqlf0TFuMWqKALnYf95bbtJX8beAerHflwCU9WJp6Sb5SvIH6pc9dmZtndqo/6Gzeq93/WzyFqpmvmAvW+1jyk6pBlvT70+HHA6SsgwJRkdgF8/juozut/bv+LKvZ62alTdXyi09Zp4Mkv82jvmiB19uMWqizc0n2qkulcC0X5HGXxraq0tKWoKTMnd6gcjZM7Bics+oQoM1D2e/btao0GFT1ljjKyhJTKQZ105sg1NLPWMOMy5wrVvoQmKp9h6T7lFzrjAdfM05eldymh/tHPVZVbW+izVEFHVwUR2MHpLSCEgIv+pUxEm77f2+CjuxP2PKPUY1fXLNKwm/cyS/nrx8fYdWIYtXTGmjN/pjYha347PGetvaRdokpQWHJWGw3w5f/ByxcrR/ztW1X9qIGmuJnroaFo6P4d0Fte++1brGst1bnKST8S85KZOdeoBLsF3x35PezhzJ+pHuuXPW1f2RNH0elUx0Zjl8rUtiSUpVTPn+fOUSXxn10NlUdcvzYLnN4CAiAwBtb+TSX07H5SHfv2HWgudzwcUcOpfHhIOffKG+wodTJecPeC69+ARTe76P6esOAmFQ7at2BgU4WpQOMjKlLrji+tR+akrjWZmd61fL4vGU+oTn5FGSrnxxLm8hojcVCbCZoEd++0XB7cmUTPVpFvkdNdO09fwqaqMvh5n8PBDf3PtdXDm9+Bj36quhtevUGZu585W4VJj7LvQhMQoPwMqetg6x9U5ETG46pw2NRBfY00xoia5o4ezaGicYg6OqcbC7+nfu5/Sf08sa1/gcbLn7GdKe4TDMlrhjYzFe9VtavOeVgljW61YJo1zx+S6PwKxqcSi29TGdYf/1Jl34Mq4vifdBWGf96f4Lr/KX/RXTtUK9PN98G7d6jWraOEJiDAZGr6pzIBbFivwvOW36NFJo0jPs4uxyjBXSc0ATGQ4Hhlxtr/strkbFgPvqFwxzb7CzTOXK/CcEv2WR+z6zHlu5h3A1z8L5Vk+N7d/XsvG7qUryPpLEff1amN2dQkjaqr4K7H4YXzQKKqMKy4t/fvFhAF39kIZ/8Kvn1bdTXUHxqdZY7KLBOBgChY+3fVVtAvQmkVGuOGDw/pSYrwIzU6QBMQllh8qyp7vv3vMP8GuH2bKqFtL6kXqqxna9FMtQWqlMmiW5Q/IyBafV9K9vbmEoDK7O5scsy8dLoQMkUl2+Z/qVrnTrtAVQa21IFO5wZnPgg3va+c1s+do1rWuhitmmtfZl2hvgiRacNrDqLhUqqaOtidX8O9ZyeTXdZIWYMmIAaRtFpFyMQtgjkj2Nx4B6nijjnvwXl/HBw5tOdpVRl4yZ29x2ZfqcZ/8UeVMBqRqsxLQqdqX2kMzaJbVImekCkqr2UobW/KSmVy2nina4IeBqAJiL4IAWf+dKxXoTEAs3lp3ZxYqls6ySyuH+sljT90Orjw/xy7x8z1cGyLyj9I6JN70FavdquzrlBBHWbMptknlipT0y2fKgd17HwVPqsxNEIogTwc/MLhhrdHxQSumZg0xj0fHiojOdKfaVH+RAV4U9PSSUf3BMmFmEhMu0D1xhhoZtr/EnS1WI7q84+EdX9XpqVtf1Q+DEfCWzXsY5T8o5qA0BjXVDa1s6eglnWzYxBCEBWoquRWNU2gUNeJgneg6o2R815vjoOhS/WwSDzTepjszMt7e2ZLg+agPoXQBMQYIKXkO8/v4Z39w+xxPU6paGzngn99zcGiOqff++Nvy5ES1s1Rpo2oIO+eOTVcwMz1qkdG8R71/+yNKultuY2cICFg3T9Ua1sPXxWSqXFKoAmIMSCvspntudVsyrJRrmAC8eLOkxwtb2KzC97PB4f0TIvyZ1qUiuOPCjALCE2DcAnTzlfhq+YKr7seU0UHk8+xfZ1fuGrPedmTWi+UUwjNST0GbM9VCV8Hi+owGiU63cTNt2jp6Oa1PYUA7Mh1bgmMisZ2vjlZywNrpvUci9Y0CNfiFWCKZtoE09eqpkQXP2pfPSRNczjl0DSIMWBHnnqQNrV3c6KqeYxX4xhv7Sumsb2bC2ZGk1vZTLkTQ1A/Oqw3mZd6q+mG+Hrg6aajXBMQrmPmelVqZvN9ymw055qxXpHGGKEJiFGms9vI7vwaViaHA7C/0Pl2+9HCYJS8sPMkCxKCuW9NMtAr/JzBh4f1TI8OIDmyt0yEEILIQC8qNROT65h2gTIz1RepXiinQE5QZ7cRg3Fi9GAYT2gCYpQ5WFRHa6eBG5clEOLrwQEXOHZHi89yyimqbeW29CTSogMJ8/Nkp5MERHlDO9+crGPd7JhB56ICvZ2qqWgMwMu/1xex+LaxXo3DdHYbueKpXTzwRuZYL2XCofkgRpmdedXoBCyfGs78hBAOFE3cpK9ntxcQH+rD+TOj0ekEZySHsyOvGiklwsE47S2HVeXWtXMGC4joQG+OlDc6dH+NIbjwb6o/gl/4WK/EYR77IpfDpQ3UtnSO9VImHJoGMcpsz6tmbnwwQT4eLJwcQl5lM/WtE++De6Cojv2FddxyRiJuJif7yuRwqpo6OFbheLXJDw/rSYsJZGqE/6BzkYFeVGgahGsJiFId5yY4h0rqefLLEwR4uVNa30ZzR/dYL2lCoQmIUaShrYus4nrSTf6H+QnBABycgKUjnt9eQIC3O1ctiu85tjJFvS9Ho5nK6tvYX1jHRRa0B1AaREunQfuya9iko9vAT97KItzfk4cuUa1Ec52weTmd0ATEKJJxogajhJUpEQDMnRSMTsDBceaoPlHVbLOtZ3FtKx99q+f6pQn4e/VaKWODfUiK8OsJ4x0pPeYlC/4HUD4IQPNDnKbUNHdQaUcU278/z+V4RTN/uWIOiyar2lC5FRM7anC00QTEKLIjrwpfTzfmxSvNwc/LnbSYQPaPI0d1eUM75/3zay57YqfVENwXd55EJwTfWzFl0Ln05HD2FtQ6VCvp8yMVTI8OIDHcz+L5SFO5DXseEhqnHg++fYhVf9vGm/uKkVYaHGUW1/P0Vye4ZlE8Z6dGEh/qi5e7juOaBjEsNAExiuzMq2FZUhie7r2/9gUJIWQW1Y+bELzssgYMRklBdQsXP7aDjQf7lwNpaOvijW+KuGhODDFBPoOuX5kSQVuXgQOFIzObdRuMZBU3sCwpzOqYaLMGoQmI05Lcyma6DZIH3z7Ej97MomWAqbG9y8CP38wkOtCbX12kemK46QTJkf4cr9Q0iOGgCYhRoqSulYLqlp78BzMLJgfT0mkYNzubo+VqHR/en86s2CB++EYWP30ri7ZOpRG88U0RLZ0GbktPsnj90qRQ3HSCHXlVI56/rcvAwsnWy0WbTUxauY3TD4NRUlbfxq0rE/nRudPYlFnKxY/tIKesN6rtn58d50RVC/935RwCvT16jk+LCnCZD6K2pZOXdhZw5VO72JBx0iVzjAWagBglzI7b9JT+AmJhQigwfhLmjugbiQ/1ITnSn9duX8p9q5N5+0AJlzyuvoQv7jzJsqRQZsUFWbw+0NuDefHBI3ZUm38PC2wICD8vdwK83LVyG6chlU3tdBslk0J9uX9NCq/etozmjm4ue3Inr+4pZH9hLc9sz+f6pQmkm3x9ZlKi/NE3tNPY3uWUtXQZjHyaXc6dG/ax9JHPefj9HA6VNvDizpNOuf94QBMQo8T2vGqiAr1Ijuwfthkf6kO4v+e4SZg7Wt7E9OhAANzddPz4vFT+e8sS6lo7ueix7egb2rndivZgZmVyOIdKG0YUvnugqI6oQC9ig2xn70YGerlMQLR0dPPO/pJxY/bT6KW0rg2ASSHKvLl8ahhbfpDO0sRQfrXxW258bi+xQT78cu3gdqvTTBn5jjqqq5o6+P37OSx7ZCt3bNjP/sI6vrt8ClvuT+fX69LIr24hf4KX0DGjCYhRwGiU7Mqr5ozk8EEJZEII5ieEcHAcJMy1dxnIr2omLTqg3/H0lAi2/CCdlSkRLEgI5uzUSJv3SU8JR0oVtTVcDhTVsSAhZMhEu+ggb5cJiKe/OsGP38rinQOnRjn2U4nSepOACO71f4X7e/HyzUt48IJUvD10/O2qOf2i68yYKwI7YmbalVfN2ke3s2H3SZYmhfL8TYvI+MUafnPRDGbEBrJ6uvpufHG0csRzjCdcKiCEEMFCiLeFEEeFEEeEEMuFEA8LIUqFEJmm19o+438hhMgTQhwTQpzvyrWNJjn6RupauwaZl8wsnBxCQXULNc1ja1PPq2zGKGF6TOCgc5EB3vz3liW8c/eKIavPzo0Pxt/Lne3DLLtR2dROcW2bTf+DmagAb5f4INo6DbyyW1Wn/ednx22G+2qMPiUmDSIupH+AhE4n+P5ZyRz4zbmsmGr5ezYpxAcfDzeOj0CDMBgl//jsODc8v4dAb3fev28lT96wkDVpUXi49T5GJ4X4Mj06gM+PVAx7jvGIqzWIfwMfSymnA3OBI6bj/5RSzjO9tgAIIWYA1wIzgQuAJ4UQbi5e36hgzgs4I9nyB3dBgnogjrUWcUSvHH3TB2gQfbGnhIaHm45lSaHD9kOYI5/mJ9ghIIK8qWxqx+hkM9DbB0qoa+3ix+dOQ9/Qzn8zTjr1/hqOUVLXRoivB76elqsE2fp86kyRTLmVw9MgKhrbueG53Ty6NZf18+PYfO/KHjOsJVZPj+Sbk3U0tDnH1zGWuExACCETRgKCAAAgAElEQVQCgVXA8wBSyk4ppa0n4KXA/6SUHVLKAiAPOCUKzO/IqyI1KoDIAMt29TmTgnDXiTH3Qxwtb8LbQ8fkMMv5B8NhZXI4RbWtFNW02n3NwaI6PN10zIqz/uUzExXgRZdBUuvEMiVGo+SFHQXMnRTEvauTOSs1gie2nTglvuinCqX1bUwK8R3x9SlR/sOKGPz6eBVr/72drOIG/n7VXP5x9Tz8LJiv+rImLRKDUfL18ZFF8o0nXKlBJAFVwItCiINCiOeEEOYnz71CiENCiBeEEObtYhxQ3Of6EtOxCU17l4FvTtb1lKGwhLeHGzNjA8c8kulYeRPTogJ6ais5gjlbfDjlv/cX1jErLhAv96EVR1c0Dtp6tJKC6hZuTU9CCMGD50+nsb2Lp7864bQ5NPqz7WjloFwbW5TWtRIXPDj/xl6mRQVQ0dhhl9B/9ut8vvvCXsL9vdh87xlcuXCSXXPMiw8h1M+TraeAmcmVAsIdWAA8JaWcD7QAPweeAqYC8wA98P9M4y09lQbZD4QQdwgh9gkh9lVVjX8Jvbegls5uo00BAcqscqikgW6DcZRWNpij5Y02zUvDYWqEHzFB3nbnQ3R2GzlU2tBjbhuKyEDnC4hnt+cTF+zD2lmqQdGM2EAumxfHCzsKtLIeLkBKycPvZ/PXj4/ZPb60vm2Q/2E4TItSUYRDOaq7DEb+vTWXVdMieO+eM0iJsv974aYTnJUawZfHq8b0++wMXCkgSoASKaWp+zlvAwuklBVSSoOU0gg8S68ZqQSI73P9JGBQk2Mp5TNSykVSykUREREDT487duZV4+EmWJoYanPcgskhtHUZehLVRpuqpg6qmztt2laHgxCClcnh7MyrsStcNEffSGe30Wb+Q1+cnSx3qKSevQW13HzGFNz7OB1/dO40jFLy763HnTKPRi/ZZY0U1rSib2intXPowos1LZ20dxl7QlxHQoop1HUoR3VWcT3NHd1ctzgeH8/hu0LPSYuivrVrQpfzBxcKCCllOVAshEg1HVoD5Agh+lZgWw98a/r3ZuBaIYSXECIRSAH2ump9o8X23GoWJIRYdaqZMUfujJWZ6aipv8L0GOdoEKCquza0dfFtacOQY3sS5OzVIAJUPSZn7eyf216Av5c7Vy+O73c8PtSXG5dN5o1visnTyjQ4lQ9NRRkBCqpbhhxvzoFwxMQUF+yDr6fbkH6IHXnVCKHyLEZCeko4Hm6CrUcntpnJ1VFM9wGvCiEOoUxKjwB/FUIcNh07G/ghgJQyG3gTyAE+Bu6RUk7oGMPq5g5y9I1Ww1v7EhvkTVSg15g5qo/q1RfGWRoE9EZt2eOHOFBUR1ywT49vYSg83HSE+3tS2eS4gCitb+PDw3quXRzfrzSDmXvPTsbX052/fXLU4bk0FFJKPjyk73nY2yUg6i2HuA4HnU6QYkck047caubEBRHs6zmieQK8PViaGMYXRyZ2PoRLBYSUMtNkDpojpbxMSlknpfyOlHK26dglUkp9n/F/klJOlVKmSik/cuXahssRfSM3vbCXwyVD74ZBfQHM8fQrU4Y2hQkhWJAQMmYC4kh5I5EBXoT6jewLYYlwfy/SYgL5yo5ojoOFdT39MewlKtA5uRAv7zoJwPfOmGLxfJi/F3esSuKT7IoxjzQ7Vfi2tJGi2lbuWKWy8guqhhYQJXUqIm5S8MijmABSogJsmpia2rs4WFw/pN9wKFZPjyS3snlYkXzjDS2T2k5e3VPIV8eruOKpXby0s8BqmWGAxvYu7n3tIP/6PJc10yOZbaVu0UAWTg6huLbNKbvi4XJU32QxQc5R1s6KZm9BrU2VXt/QRllDu93mJTPO6E3d1N7F63uKuHBWtM3wyVtXJhLu78VfPjpq82+vYR8fHC7DXSe4dF4ssUHedpuYArzcCfRxrFPytCh/qpo6rJaC2Z1fi8EoWZnsmI9zTZrKqp7IZiZNQNiBlJIvjlSyMjmcVdPCefj9HO7csJ+G1sGhcodK6rno0R18nF3Ogxek8ux3F9kdNmpOEBtpqeyR0m0wklc5uMSGM7hh2WS8PXQ8v73A6hjz+7Ung7ovUYHeDgvTN74ppqmj22p1WjN+Xu78YE0yewtq2XZsYpsNxhqzeWllSjjBvp4kRviRb6eJKS7Ex+F+5+aIJGtaxM68anw83FgweXga7UAmh/mRHOnP1glsZtIEhB0c0TdR1tDOJXNjefa7i/j1ujS+OFrJ2ke395gcpFRJVlc8tYsug5E37ljG989KHrIsRV9mxQXi6abj4CibMQqqW+g0GJ3qoDYT6ufJFQsmsTGzlKomy+agA0V1eLnrSBumBhMV6EV1cyed3SMLJew2GHlx50kWTwnpaeJki2uXJBAf6sPTX+WPaD4NxaGSBkrq2lhn6hiYGO5HflXzkJpZSV2bQw5qM9N6BIRlrXZ7bhVLEkPtyscZijXTI9lTUEOTkyrIjjaagLADc8LL2dMjEUJwW3oSb9+9AiHg6qczeGJbHndu2M/vP8hhVUoEW+5PZ9EU22GtlvByd2NWnEqYMxqlxZcrOFLufAd1X25ZmUhnt5ENJp/MQPYX1jFnUlC/Rkr2YG4cVDXCGlYfZ5dTWt82pPZgxsNNx7WLE9hbUGuXSUTDMh8e1uPhJjhvhso3SQz3p7G9m9oW21nxpXWO5UCYiQ3yxt/L3WIuhL6hjRNVLXYFltjDmrQougzS4Ta8Y4UmIOxg69FK5sYHE2EKrQSYFx/Mh/enc05aFH/75BhfHK3k1+vSeO6mRYQ44OhdODmEfYV1JP1yy6DX1F9t4UdvZg7qoOUoR/WNuOsEUyP8hx48AqZG+HNOWiSv7C4cVPyuvctAdlmD3fkPfXGkN/Xxiib++vExpoT5ck5alN3XXblwEjoBb+8vHnrwEHyaXc6KP2/lk+xyh+81UTCbl9JTIgjyVRFjSRGqwIItodvQ1kVTR7dDORBmhDB1l7NgYjI/yB11UJtZkBBMkI/HhDUzOebtOQ2oauogq6SeH50zbdC5IB8PnrpxAR8c0pMY7me1ic5wuHVlEkE+HlhKwKxqbue1PUVkFtfzxPULhm2SscbR8iaSI/2HvYMfDrelJ3HtM7t590Ap1y9N6DmeXdZAl0EO20ENvQJiOL2ppZS8ua+YhzZn4++l/n7DKS0SFejNWamRvL2/hB+dm+pQWZL3Mkspa2jnzg37+d6KKfxi7XSnmDXGM5nF9ZTWt/Gjc3u/T0mm3uP51S1WNe/eHAjHIpjMTIvyt1iSe2deNeH+XqQOI3PaFu5uOpVVfawSg1E6pYzNaKIJiCHYdqwSKWF1muUeCEIILp4b67T5ooO8uXd1itXz62bH8oP/HeTSJ3by0MUzuH5JgsNOu6P6RhYPkentKEsTQ5kVF8hzO/K5dnF8j29muAlyfYkKNCXL2Skgmju6+fXGw7yXWcaKqWH869p5Vgso2uLqRZO465VKvj5exdnTbffGsIbBKNmZV8Nl82IJ8fPkxZ0n2V9Yx+PXz3dKscTxyoeH9Hi66ThnRq/WFhfsg4ebsKlBOCMHoi/TogJ4c18JtS2dPaHdRqNkZ141Ky30bXGENWlRbMosI7O4ftiBGGONZmIagi+OVBIT5M0MF4SAjoSBHbTue/2gQw6whtYuyhraXeZ/MCOE4Pb0JPKrWvjyeO/O7UBhPfGhPv3Md/YS4uuJh5uwKxcip6yRSx7bweasMn507jQ23Lp0RMIBYPX0KML8PHlz38jNTN+WNtDQ1sXZ0yN56OKZ/Oc7CymsaWHdozv44NCgCjOnBEajZMthPaumhRPk05uQ6O6mIyHU12YXNnMOhDOc1NA3kqnXD3G0vInq5k678paGw5kpEbjpBF9MwHDX01JAGI2S3flDdzvr6DawPbeK1Sbn9Hihbwetj74t56LHdthVzsISx0xfEFdEMA1k7ewYYoK8efZrFfIqpWR/UR0LR6A9gMqKjQyw3VnOnLB42ZM7aens5rXbl3H/mhSHVH1Pdx3r58fx+ZGKETd5MmeXm7PNz58ZzZYfpJMS5c+9rx3kVxsPn3LNig4W11PW0M66OTGDziWG+9vWIOra8HJX2fPOwFLRPnNhyZVW+raMlCBfDxZPCZmQfojTUkC8ua+Ya5/ZzTcna22O25NfS0unoSfhZTxh7qD1vzuW0dlt5LpndlvMyxgKcw2mNBdrEKCigL63YgoZ+TV8W6pCHauaOkbkoDYTNURv6i+OVvLr975leVIYW+5PZ1nSyGrrDOSqRfF0GSTvZY5st78jt5q0mEDC/Xs1p0khvrx553LuXJXEq3uK+MtHp1Zpjy2HlXlpjYWggKkRfpysabVa2NFZORBmogO9CfBy7+eo3p5bTUqkv93lXobDmulRHC1vmnD1vE5LAXHpvDiiAofOit16pAJvD53VFobjgcVTQnn+psU0dXTz2t6iYV9/RN9EsK9Hjz3f1Vy7JAE/Tzee31HQk0MyEv+DmaF6U7+XWUaIrwfP37SIMH/nvcfU6ADmxgfz1r7iYWdWt3Ua2F9YZzGU0sNNxy/WpnHVwkn875siq9m+E41e81KExXpXieF+dHYbKTP5GgZSWu+cHAgzQoh+zYNU35Zaq10fHWX9gjg83XW8sNN6wuh45LQUED6ebjxwzjT2F9bxWY5lu6CUkq1HVfa0t8f4jiyZERvIyuRwXtpVMOykMXMPiNEyoQX5eHD14njezypjy2E9Ph5uDvWgiLTRm7q9y8DWIxVcMCumXwlvZ3H1okkcLW/ikJ31uczsKaih02C0acq4LT2J9i4jr+4ZvtAfjxwsrkPf0M5FFsxLoAQEYDWjuqSuzSkhrn2ZFhVArmlHf6CwjvYuo9PyHwYS7u/F5fPjeGd/yZj3nh8Op6WAALhq4SSSIvz46yfHLDb1yK1spqSujdXT7Y+RH0tuTU+korGDDw/bb/IwGiXHyptc7qAeyC1nJGKUkk+yK5gbH+TQwzs6yJvmjm6aLeSGfHmsktZOg9WHkqNcPDcWbw/dsJ3VO3Kr8XTXscRG5FhqdADpKeG8tOskHd0T3xfxwSE9nu46q+baRHMuhAVHdWunSqJzpNWoJVKiAqht6aS6uYPtedW46wRLnWSCtMStKxPp6J5YQv+0FRDubjoePD+VvMpm3j1QOuj856bs6dUjDGMcbc6aFkFKpD/Pfm27kGBfiutaae00OK2LnL3Eh/py4Sz10HbEvAS9oa6WzEwfHNIT5uc5ZLOmkRLo7cHaWTFsziyjrdP+h/iOvGoWTwkZUjO9PT2JqqYO3s/S2xw32nxzsnZYQstsXjprWgQBFsxLABH+Xvh7uVt0VJvNTs40MQGkRCpH9fGKJnaY+rb4D9Fv2qH5ogI4KzWC/2actCsAoaWjm4wTNWNaHPK0FRCgIkfmxQfzj8+OD/qDfXGkkllxgS5xWLkCVQIkkRx9Ixknho7QAuV/AFxSxXUo7liVhLtOkO5gSGGUldajbZ0Gth6p5IJZ0S4xL5m5alE8TR3ddmdDVzV1cLS8yS5bd3pKOKlRATy3PX/cVJAtrGnhqqczeOMb+7WmHH0jFY0dXGBq5WoJIQRJVor2ldQ5NwfCjLkm0zcFdXxb1uAy/0Nfbk9Porq5k02ZgzelfZFS8qM3M7nu2d3c+9pBGseoltNpLSCEEPz8wumUN7b39AQAqG3p5EBRHWsmiHnJzKXz4gj39+S5HfY5wo6VNyFEb8jfaDI3Ppish84bcccuM9YExLZjlbR1GSyGVDqTpYmhJIT62v3A3GkKb023o5S0EIJb0xM5Wt7Ezjz7hL6rMftbsort97tkl6mx84fQFhPD/SxqEGYB4WwfRFSgFwHe7ry6pxApnVdewxYrpoaRFhPIc9tta/qbs8r4JLuC9JRwPs4u56JHd3CoZPTbl57WAgJgWVIYZ6VG8MS2vJ4w0S+PVWKUjMvwVlt4e7jxnWVT+OJoJXlDdMwC5aCeEuY3ZDtUV+HnBHXeWm/qDw/pCff3ZGmi62zKoMKNr1o4iYz8Grsaw2zPrSbE14OZsfZpbZfOiyUiwItnt4+PCrLZZY2mn8MREI34e7kzOdS2DyEx3I/S+rZB2nxpfRvuppwXZyKEYFpUAJVNHQR4uzN3kuOlcuyZ87aVieRWNlttpFXZ2M5vN2UzPyGYl25ewht3LKPLYOSKp3bx4hC9aJzNaS8gAB48fzpNHd08+VUeAFuPVBIR4MWsWNd/YJzNjcsS8HLX8fyOk0OOPVreNOr+B2fj7+WOv5d7v4J9rZ3dbD1awYWzYkal9s2ViyYh7CjgJ6Uq5bAiOdzuMvBe7m7ctHwyXx2vGrKP8mhgFgx5lc12+yGyyxpJiwkY8j0nhvshJRQOELSldW3EBHu75G9p1p6XJ4W51BTZl4vnxhIV6MVzFnqkSCn5pSlJ8u9XzcVNJ1g0JZQt96ezKiWC39noReMKNAGBChNdPy+Ol3aepLi2la+PV7FmeuSwejmMF8L8vbh8wSTePWA7nK61s5uTNS2kTnABARAZ6NWvcdAXRytp7zK63LxkJibIh1UpEby1v8RiRJyZE1XNlDe2DztT94alQzddGglSSl7bU2Qzj2Tg+JyyRkJ8Peg2So6XD530ZTBKjugbmWnHZispXD2sC6r739fZORB9SYlUn//RMC+Z8XTXcdOKKezIq+aIvrHfuY0HS/n8SCU/PT+1X3XlED9PnrvJci8aV6IJCBM/PHcaUsLt/91HU0f3hIlesoQ5nM5a/wVQ3bSkdF0PiNEkekDr0Q8P6YkI8GLxCHpyjJQbl01G39Bu0//TU0p6mAIixM+TKxdOYuNB602XRsLOvBp+ufEwL/Xxv9miorGDmpZOLpsfB9hnZjpZ00Jrp4EZdpjUzKGuAx3VJXWtTg9xNZOeEs7UCD/OnTG6/sbrlyTg4+HWT4sob2jn4c3ZLJocws1nJA66ZmAvmq/t6PXuKJqAMBEf6suNyyZztLwJT3fdqO4onE1ypD+rp0eyIWNw/wUzR007l7RRqMHkaqICe5PlWjq61Q5rVvSollY+Jy2SC2ZG849Pj1tsRAMq/2FKmC/xQ9jiLXHLGYl0GY1syDjp2EL78NwO5dc4UGjfTtQsEC6cFYO/l3uPP8L2NWqMPT4Xfy93IgO8KKjqFRCd3UYqmzpcp0FEBbD1x2cRE+Sa+1sj2NeTqxdNYnNWKZWN7Ugp+cW7h+g0GPmbybRkjXnxwWz5QTr32aj67Cw0AdGHe1cn4+/lzhlTw8bMcessbktPpKalk/cO9obTSSn5trSB372fzf99fJRAb3fiXbQzG03MvanN2e8d3UbWzXFeCXZ7EELwx/Wz8Pd258dvZQ0yNXUZjOzOrxnxxiMpwp8106PYYKHp0kjIrWjiy2NVBHi5k1VST5cN05iZ7LJGhFAm2RkxgXZpENllDXi4iR5TzlAkhvcPddU3tCGl80NcxwM3n5FIt1HycsZJ3tpfwrZjVfzsguk9WeW2CPT2GJUNkCYg+hDq58k7d6/gz5fPGeulOMzypDBmxATy3I4CKhvbefbrfC7893YuemwHr+4uYvnUMF66ZcmE9LMMJCrQiy6DpLalkw8PlREZ4MWiMai7H+7vxR8uncWhkgb+83X/qKPM4npaOg0OVQq9PT2RutYu3jlQ4uhSeW57Ad4eOn56QSrtXUaO6od2gGeXNTAlzA9/L3dmxAZyRN9ktbhezzWljaRGB9jdjCopon+oq7lR0CQXaRBjyZRwP86bEcUru4v4w/s5LEkM5ablU8Z6Wf3QBMQAUqMDJkxynC2EENy+KpG8ymaW/nkrf9pyBC8PN/5w6Uz2/moNT96w0OEs5vGCuTd1fnULXx6rYu3smDETfOvmxLBuTgz/+vx4T6VcUP4HnYDlDhR+XJIYypxJQTy/vcCuHb81qpo62JhZyhULJvVUVrXH4Zld1tjjS5gZG0hbl2GQQ7kvUkqyyxqYGWN/NGBiuB+1LZ09RQp7cyAmvqZridvTk2ho68IgJX+/cu6427BpAuIUZt3sWC6bF8tdZ07l8x+tYtM9Z/Cd5VMI9nVOTf3xQqRJQLy2p8hkXhqd6CVr/OHSWQT5ePDjN7N6HuQ7cquYMym4X6Oc4SKE4L7VKeRXt/D0lydGfJ8Nuwvp7DZyy8pEYoO8iQr06unsZ4361k5K6tp6fAnmqCRbfgh9Qzt1rV3MjLM/EKI3kklpESX1bQjBKbFps8TCySHcuSqJf1w9l4Sw8ScENQFxCuPpruNf187nZxdMJ9lOG/BExPzw+OBQGVGBXiNuQOQsQv08+eNls8kua+TJbSdobO8iq6TBKZVCz50RxSVzY3n0i1xy7HASD6S9y8Aruws5Jy2SqRH+CCFYkBAypAZhnsucG5QS5Y+nm86mgBiOg9pMT9E+k4AorWsjKsDbpf3SxxIhBL9Ym8YFs8Z2U2MNu37rQoipQggv07/PEkLcL4QIdu3SNDTsI8LU56HLIMfUvNSXC2ZFc+m8WB77IpcXdhRgMEqn1fr53SUzCfLx5CdvZQ27vPu7B0qpbenktvSknmMLJ4dQUtdGpY18iIEPew83HanRATYd1dllDQgxvFDq+BBf3HSCfFMkU0ld6ynpoJ4o2CuW3wEMQohk4HkgEXjNZavS0BgGnn1aUbqqtPdIePjimYT4efKvz3Px9XRzms8nxM+TR9bPIkffyBPb8uy+zmiUPLcjn1lxgf0q3JprJNnSIrLLGogO9O7XdGlmbCDZZY1WSz9klzWSGO43rJIqnu464kN8ejWIeuf3gdCwH3sFhFFK2Q2sB/4lpfwhMH6+iRqnPZEB3sQEeTM/fvw43tWDfDagivo500xy3sxo1s+P44lteXb3I//yeCX5VS3cnp7Ur0HUrLhAPN10HCiyXgwuu6xxkKloZmwg9a1dlDVY1jxyyuzLoB6IOdTVYJSUN7S7LAdCY2js/cR2CSGuA24CPjAdG7m3TUPDyfz0glT+csWccWFe6su5M6J4ZP1sHjhnmtPv/dDFMwj1U6Yme+oiPft1ATFB3qyd3X9v5+Xuxqy4QKsJc22dBk5UNQ8SEDPMjmoLAqqupZPS+jZmDcP/YCYpwp+T1S3oG9roNkrNxDSG2CsgbgaWA3+SUhYIIRKBV1y3LA2NwXQbu3kv7z26jYO7x52dGsmZ0xzrLeEqrl+awNx457vsgn09+fPlszla3sRjW22bmr4tbSAjv4bvrZiCh4WidAsSQjhU2mDRp3G0vBGj7BUIZtJiAhDCciRTjt7ssxiZBtHWZeiJrNI0iLHDXgGRBDwgpXwdQEpZIKX8i+uWpaExmK+Kv+I3O3/DXv3esV7KuGFNWhRXLpzEU1+dIKvYuono+R0F+Hm6ce2SBIvnF0wOobPbaNHpbC0aydfTnaRwP4vXmM1ew4lgMpNkyiQ21646VXMgJgL2CohrgVwhxF+FEGmuXJCGhjUyqzIBqGyrHOOVjC9+c9EMIvy9+MlbWVQ2tlPT3NHvlVvRxPtZZVy9ON5qHsbCyWZH9WAhk13WSJCPh0Vn8czYIIsaRHZZI7FB3oT4DT/nxhzqusMkIDQNYuywK7xASnmjECIQuA54UQghgReB16WUY1+kXuO0ILNSCYiatvHRXW28EOTjwV+umM33XvyGJY9stThGJ1TBP2tEBXoTF+zDgaI6bqX/uJyyBmbEBPZzbJuZGRvI5qwyals6Ce0jDLLLGgaZpOwlKsAbHw83yhvbCfPzxMfTdu9uDddhd/yZlLJRCPEO4AM8gIpo+qkQ4lEp5WOuWuCpys++/hnnTD6HcyefO9ZLmRB0GbrIqckBoLqteoxXM/44KzWSDbcusdiyE2BymN+QVWTnJwQPyqjuMhg5Ut7ETcsnW7xmVpw5o7qhp794a2c3+dUtXDTCgok6nSAx3I8cfaMW4jrG2CUghBAXA7cAU4ENwBIpZaUQwhc4AmgCYhg0dDSwpWAL3cZuTUDYyZHaI3QaVX0eTYOwTHpKRM9DeiQsnBzCB4f06Bvaespfn6hqprPbaNXZbPYxZJc19sx9RN+ElCPzP5hJjFACQotgGlvs9UFcBfxTSjlHSvk3KWUlgJSyFSU4NIZBQYNqEpJbnzvGK5k4mM1Lk/wnUd2uaRCuwJzId6Cw1w+RXWq7XEawrydxwT79/BA5Jqe1WbsYCWZH9Xj0P7R1t7ElfwtGOfKCiRMFuwSElPK7UsqvrZyzbPTUsMqJelVoraixiA6D8zqEjSUnG066tJl6VlUWMX4xpIWlaRqEi0iLCcTLXdfPzJRd1oi3h46kPu0vBzIjtn9viGxTW9IYBwrsJY5jAbEpbxM/2/4zvir+aqyX4nLsrcW0TAjxjRCiWQjRKYQwCCGGXylMA4ATDUpAGKSBkw0nx3YxTuBE/Qkuee8SPsj/YOjBIySrKot5EfMI8w7TfBAuwtNdx9xJwf1KbmSXNTA9OtBmc5qZsYEUVLfQ0tFtukZlUFtyatuL2aQ1bRz2TDdH071x/I0xXonrsdfE9DgqgikX5aS+Dc3vMGLyG/Lx91A7suN1x8d4NY6zo3QHEsn2ku0uuX95SzkVrRXMjZxLuE84jZ2NdBo6XTLX6c78ycFklzXQ3mVASkmOfnCJjYHMjA1CSpVQ12Uwcqy8ySH/A6i+LNsfPJvlSWEO3ccVZFVm4Sbc2Fm6k+LG4rFejkuxuziMlDIPcJNSGqSULwJnD3WNECJYCPG2EOKoEOKIEGK5ECJUCPGZECLX9DPENFYIIR4VQuQJIQ4JIRaM/G2Nb/Lr81kZtxJ3nTt59fYXWxuvZOgzANhTvscldlnzjm1uxFzCfNQDo7a91unzaCg/RJdBtaYtrm2jqb17yGxoszD4trSR3IpmOg3GnsZCjhAf6uuQFuIKqtuqKWku4fq063ETbrx5/M2xXpJLsVdAtAohPIFMU7LcD4GhG6fCv4GPpZTTgbmoiKefA1ullEAHnwcAACAASURBVCnAVtP/AS4EUkyvO4Cn7H8bE4eWrhb0LXqmhUwjKShpwguITkMn+8v3E+ETQW17rUs0oqzKLLzdvEkNTSXcR5XM1sxMrmFBn8quZr/CUNpATJA3Ib4eZJc19Llm5A7q8UxWVRYA500+j9UJq9mYt5H2butl0ic69gqI75jG3gu0APHAFbYuMCXWrUKVB0dK2SmlrAcuBV42DXsZuMz070uB/0rFbiBYCHHKVYw1RzAlBSeRHJxMbt3EjmTKrMyk3dDOXXPvAiCjLMPpcxyqOsSMsBl46Dw0AeFiIgK8SAj15UBhPdlljbjpBKlD+AGEEMyKUxnV2WWN+Hi49TiZxwvZNdnk1Tm+GcuqysJd505aWBrXpl5LQ0cDnxZ+6oQVjk/sjWIqBAIALynl76SUPzKZnGyRBFShMq8PCiGeE0L4AVFSSr3pvnog0jQ+Duhr0CsxHeuHEOIOIcQ+IcS+qqoqe5Y/rshvUM3sk4KSSAlJQd+ip7nTel/f8U6GPgM34cbaxLUkByc7XUB0GDrIqc1hbuRcAE1AjAILEoLZb9IgUiL98fYYOpN5RmwgxyuayCqpZ0asbaf2aKNv1nPrJ7fy571/dvheWZVZzAidgZebF4ujF5MYlMgbR09dZ7VNAWHyCzwshKgGjgLHhRBVQojf2nFvd2AB8JSUcj5K8/i5jfGWPlGD4iallM9IKRdJKRdFRIzP6p22OFF/AnedO/EB8aQEpwBMaDPT7rLdzImYg7+nP8tilnGg8oBTQ3dzanLoNnYzN0IJiFBv1ehGC3V1HQsmh1DV1EFGfo3dvoSZsUF0GSQHi+oddlA7EyklD+16iJauFkqaShy6V5ehi+yabOZEzAGU5nRN6jUcqj5Edk22M5Y77hhKg3gAOANYLKUMk1KGAEuBM0x+CFuUACVSyj2m/7+NEhgVZtOR6Wdln/Hxfa6fBJTZ/U4mCPn1+UwJnIK7zp3kkGRg4ibMNXQ0kF2TzbKYZQAsj11Oh6GDg5UHnTZHVqWy+ZoFhKebJ4GegZoG4ULMfoj2LusZ1APpKxTGk4B4O/dtMvQZxAfEU9FaYbFUvL0cqztGh6GDeZHzeo5dPPVifNx9ePPYqemsHkpAfBe4TkpZYD4gpcwHbjSds4qUshwoFkKkmg6tAXKAzajGQ5h+bjL9ezPwXZPWsgxoMJuiTiXyG/JJClL9gGP8YvB193WKbXQs2KPfg0SyPHY5AIuiFuGuc3eqmSmzKpNJ/pN6TEugzEw17ZoG4SqmRwfgayqQZ+/DPjHMr88148NBXdpcyt+/+TtLY5Zy86ybMUgDVa0jN0ubHdTmzQpAoGcgaxPXsiV/C42dp15q2FACwkNKOWirJqWswr6OcvcBrwohDgHzgEeAvwDnCiFygXNN/wfYAuQDecCzwPftegfjgI25GzlcdXjIce3d7ZQ0lzA1eCoAOqEjOSR5wmoQGfoM/Dz8mBU+CwBfD1/mRsx1moCQUpJVldXjfzAT5hOmmZhciLubjjmT1EPeXhOTTidIiwnEXSdIibKede0MMsoy2Ji70WZItVEaeWjnQwD8fsXvifNX7szS5tIRz5tZmUmUbxTRftH9jl+Teg3thnY2520e8b3HK0MJCFvZSENmKkkpM03+gjlSysuklHVSyhop5RopZYrpZ61prJRS3iOlnCqlnC2l3DecNzJWtHe38/vdv+eprKGjcgsbCzFKI0nBST3HUoJTyK3LdWmZClexu2w3i6MX46Hr3Sssj1nO0dqj1LVbbl85HMpayqhuq2ZexLx+x8O9wzUTk4u5elE8l8+PI9Db/s7CVyyYxDWL4/Fyd2157scPPs5vd/2W+764j/p2y02S3jz2JnvK9/DTxT8l1j+WWD9VWVbfMnKjRFZVVj/zkpm0sDTmRMzhjWNvTMjvsS2GEhBzhRCNFl5NwOzRWOB4x+xEPVh5EIPRdl/gvhFMZpKDk6nvqJ9wJpPipmJKmktYHrO83/FlscuQSPaU77Fypf0M9D+YCfPRym24mssXTOIf1wx+GNri+qUJ/Gm96x8LZS1lTA6cTEZZBle+f+Ugn1dxUzH/2P8Pzog9gytSVDR+jL+KmC9rHplbs6KlAn2LftBn0cy1qddysvGkUz734wmbAkJK6SalDLTwCpBS2r+1OIUx2yWbu5o5VnfM5tgT9SfQCR1TAqf0HEsJUZFMEy0fwmxGWha7rN/xmWEzCfAIYHfZbofnyKzKxMfdp+d3ZCbcJ5zW7lZau1odnkNjYtFh6KC6rZp1SevYsHYDnm6e3PzxzTx3+DmM0ohRGvnNzt/gLtx5eMXDPZnYXm5ehHmHjViDsOR/6Mt5U84j2Cv4lHNW211qQ8MyWVVZBHkpe+3+iv02x+Y35JMQkICnW2/nreRgFck00UJdd+t3E+UbRWJg/+5j7jp3lsQsIaMsw2F1O6sqi9nhs3HX9W9bYi63MdG0Lg3H0TerB3ysXywzw2by5kVvcu7kc/n3gX9z9+d383TW0+yv2M+DSx4c5CuI848bsQ8iqyoLT50naaGWOy57uXmxPnk9XxR9QUVLxYjmGI9oAsIBpJRkVmaSHpfOJP9J7Cu37TbJr88nMaj/AzXMJ4xQ79AJpUEYjAb26PewPHa5xVo5y2OWU9ZSRlFT0YjnaO1q5VjtMYs7NnNEk+aoPv0oa1Emolh/5VPw9/Tnr6v+ym+X/5Z95ft4KuspVk1axaVTLx10bYx/jEMaxMzwmXi4WTecXDXtKozSyLt5745ojvGIJiAcoLS5lJr2GuZGzGVR9CIOVB6wGlnRZeyisLGwJ4KpLynBKRNKgzhSe4TGzsZB/gcz5rBXR8xM2TXZGKTBpoDQ/BCnHz0ahH9vO1MhBFdNu4rX1r3GFSlX8PDyhy1uXGL9YtE364ddULLT0ElOTY5V85KZ+MB40sLSOFBxYFj3H89oAsIB+tolF0YtpL6jvqcZ0ECKG4vplt39HNRmkkOSyavPmzAdqsz+h6UxSy2ejw+IJ9YvtqfK60gw/27NWat9CfM2mZg0DeK0o7S5FJ3QEekbOehcamgqD694mAhfyxUWYvxj6DR2Dvtzk1OTQ5exa1A0nSWSg5OtPgMmIpqAcICsqqweJ+qiqEWAdT9ETwRT8GABkRKcQlt3m0Mx2qNJhj6D1JDUHl/AQIQQLI9dzl793hFnrmZVZTElcAoh3iGDzoV4hyAQWuvR0xB9i55I38h+odX2Ys6FMJup7KVnIxhpW4MAmBo8laq2Kho6GoYcOxHQBIQDZFZm9jhR4/zjiPKNYl+FZT+EeVcx0KkL9JTcmAgZ1a1drWRWZvaYkayxLHYZTV1NI6pRI6UkqzLLovYAyhEe4h2imZhOQ8qay3pyGoZLjJ8KdTWbqewlqyqLOP+4ftn81jAHnZwqWoQmIEZIa1crx+uO99glhRAsil7EvvJ9FqN3TjScIM4/Dl8P30HnzB+qiZBRfaDyAF3GLqv+BzNLo5ciECPyQxQ3FVPXUWfT5qvlQpye6Fv0/fwPw8F83XA0CPNmZSj/gxmzj3Ei+RRtoQmIEWJ2ovbNrFwYtZCa9hoKGwsHjS9oKBgUwWTGz8OPOP+4CaFB/P/2zjy6rfpO9J+fFkve18SJ4zh27CTe4mx2SIBshCZAWMuUkgmElteBKfBO+06nLT2vFB4tM5Rm3mPoMgVaSmjpkAJhoJQlEAghkJCtcTY7iRMrXuPdlndr+b0/pCtkW5IlWfKW+zlHx9K9P937u9dX93u/+/66/URoIlia6rvhX6Ixkdyk3KD8EIpK7ylrVSHFmEJrr9pV7nLCYrfQ0NPg0gQCJVofTVxEXEDJcpe6L9HY2+i3gJgZPZNIXaSqQVzuuJyoKV+aQRQ/xFAzk81uo7Kjkuz44RFMCjkJk6Mm0/76/SxJXYJRZxxx7Mq0lZQ2lQac0FbaVEq0Ptrn+UqJVMttXG409jRil3aXLyEYZsXMCkhABOJ/AEd9tez4bFVAXO6UNjqcqAnGBNeyzLhMko3JwxzVdV119Nv6PYa4KuQk5GDqMGGxWcI259HS3NvMubZzrvLeI7EybSVWu9WrX8YbX9R/weLpi9FqvNf0UQTEZK5902Pp8avIo4oD5caulM0IhpnRgeVCHGs6hlFrZH7ifL+/k52QrZqYLmdcVUaHqJ1CCJalLuNww2A/hBLB5M3EBI6SG1ZpxWQ2hWXOoeBAvcOfMJKDWmHJ9CUYtIaAqruaOkyYzCZWz1rtc1xyZDID9gG6LJO3G9/OczvZ8s4WqjurRx6s4rqxB+ukBocfoq6rzu8Hi9LGUgpTCgOKmspJyKGlr8VrIcHJhCoggqCqs8rhRPWgdhbPKOZS96VBjrDzHQ5101OIq4LLUT1BM6or2ip4tvRZkoxJXssNDMWgNbAsdRn7avf5/YP8pOYTANbMXuNznBJiO5nNTDVdNUgkH1V9NN5TmRQoYeCj1SB6rD1+haH2Wfsoby332/+gMJUc1aqACAJfhbuWpS4DGFR243z7eaZHTicuwntt/az4LLRCO+EuKiklb5x7g81/20znQCfb1mxDI/y/bNZnrMdkNvl9XHtr9pKTkDOinXkqZFMrNXt2V+0e55lMDuq76kmJTMGgNQS9jUByIU61nMIqrQELiKkU6qoKiCAobSwlRh/j0Ymak5BDvCF+kB+isqOSrATv5iVwtNLMjMucUI7qbks3P9r3I37y+U9YNH0Rr938GiUzSgLaxvqM9WiEhvdN74841jxg5mjDUdbOXjvi2KmQTd3Q4xAQxxqPTWpBN1bUdQefA6GgaB/+5EIE6qBWmBE9g2h99IR72AsGVUAEwbEmR4KcJyeqRmhYOn2pyzErpeR8+3mfETkKOYk5E8bEVN5azp1v38m7le/y4OIHefbaZ/1KFBpKcmQyJaklvG96f0Qz0+e1n2OVVtak+zYvgVvBvhBWdO239XO65TRvnHuD7ae2j9jfY7Q0dDeQn5yPRPJx9cdh3ddUoL6rflTmJfjSf+FP1YLSxlIyYjNIMiYFtA8hhCOSqWPyaxC6kYeouNNt6aaivYJriq7xOqY4tZiPqz+mobsBiaTH2uMzgklhXsI83je9T4+lx2NC3Vjx6tlXefKLJ0kwJPC7Db8LWGsYyobMDfz0wE8523aWBUkLvI7bU7OHREMiC1NGbjoTb4hHJ3QjPnnb7Db21e6jz9Y3bJ1EcqnrEmfazlDeWk5lRyU2+aVQKEwpdJkMQ43FbqGpt4lb592Kud/M7qrdfG3+18Kyr4lMc28zfdY+0mPTfY6zSzv13fWsz1g/qv3FG+KJ1EWOGMmkBKJcmXZlUPvJTsh2+dMmM6qACJATzSewS7vPwl3LZjhuKkcajrh6RfiKYFJQSm6cbz/Pwmnj07CvprOGx/c/zoqZK/j56p8H/PTkiWvnXMsTXzzB+6b3vQoIq93Kvtp9rElf4zO8VUEjNCRFJo0oID6p+YTvfPwdn2NSo1JZkLSAdbPXkZuUy4zoGWx5Zwsnmk6ETUC09LYgkcyInsH6jPW8XP4ynQOdxEbEhmV/E5UnDjxBZUcl/33rf/sc19zbjMVuCTqLWkEI4VcuRH13PS19LUH/DrMTsnmj4g1a+1pD8hsaL1QBESBKG0xfF05uYi7R+mgONxx2CQZ/NQhwlNwYLwGhtG/8l+J/CdmFnWRMYvmM5ey6uIv/ueR/eizFXNpUSkd/B6vTfYe3upNsTB7RB3Gy+SRaoeWVG19BJ4Zf7smRyR4LAs6KmcXx5uN+zyVQLnVfAhzCKSchh+2nt7Ovdh/XZ10ftn1ORM62naWqs4puSzfR+miv45Qb+mgFBPiXC6H8790TYQPB3VGdNGPyCgjVBxEgx5qOkR2f7TMiSavRsmT6Eo40HOF8+3kSDYl+3WzTY9Mxao3j6odQspiVCzxUbMzcyEXzRa9tWT+p+QSd0HFV2lV+b9OfbOry1nKy4rPITcolJzFn2MuTcADHjeFEc/iS2BQHdWpUKoumLSLZmHzZRTNZbBaXL+BMq+92vcoNPdgyG+6kxaSN6IM40XSCCE1EQAly7kyVUFdVQASAXdo53nTcZ40gheLUYi50XOBIwxG/zEvgTNNPyB7XSCalzac/Zp5AWJ+xHq3Qssu0y+P6vdV7WTZjGTERMX5vMyUyZUQNory1nPzk/IDmCg4N8VL3JZp6mgL+rj8oIa4zomegERrWZazj05pP6bf1h2V/E5HqzmqXz6estcznWOWGHgoNIi0mDfOAmW5Lt9cxJ5pPkJec57ODnC9So1KJ0cdM+lBXVUAEgMlswjxg9isuWrFdm8wmv8xLCjkJOeNWtG9ohdpQkmhM5IqZV3iMZqrurOZ8x3m/opfcSY5MpqWvxWujpebeZpp6m8hNyg14voqjPFxmpoaeBiJ1kS5NdH3GenqsPXxR/0VY9jcRqeyodL0vby33Oba+q554Q7xPM5S/KJFM3vwQFruF0y2n/QqW8IYQguyEyV+TSRUQAaD4H/y5gRYkFxCpiwT88z8ozE+cT0tfy7jExSsO+HAICHCYmao6q4bdDPbW7AVgbfragLaXEpmCTdq8ZsWWtTieSoMREHnJeeg0urDVSrrUfYnUqFSXP+aKGVcQo4+5rMxMlWaHgFg0bdGIAiIUORAKrlwIL36Ic23n6Lf1e+1H4i9TobucKiACoLSplLiIODLjM0ccq9fqXReYvyYmgIKUAsDhXB1rfLX5DAXXzL4GndANS5rbU72HrPgsZsfNDmh7I5XbUG46wQgIg9bAgsQFYfNDNPQ0kBqV6vqs1+pZlb6KPdV7wp5/MVEwdZiYHjmdZanLqGirYMA24HVsfVd9SPwPMHIuhPJQMBoNAhwPhm39bZM6mVMVEAFQ2uTocuZvqYnlM5YDBOTwzUvKQyM04yYg5sbPdYXmhpoEYwJXpA02M3UNdHG44XDA2gN8mU3tTUCUtZaRHpMedOjowpSFnGw+GZYbdkNPA6nRqYOWrc9YT2tfqyuSbKpTaa4kMz6TvKQ8rNLq1aErpXRoECHwP4DjwSJCE+E1m/p483GSjEmjKisOX1oOJrMWoQoIPzEPmKlor/CrcbnCXXl38fsNv/fYYN0bUfoo5sbP5WTL2AoIbxVqQ83GORup6arhdOtpwNFfwmq3BhTeqjBSNnV5azl5yf4VFvRE0bQieqw9rmq8ocJmt9HU0zRIgwC4etbVRGgiLgszk5QSU4eJzLhM1//Im5mpvb+dXmtvyASERmiYGTPTaz2mE80nWJiy0GM4diAoD4aTOZJJzYPwE0XtDKQuS5Q+iuUzlwe8r8KUQvZU70FK6ddFeqLpBJ/Xfc79i+4PeF8KJrOJjv6OsAuIazKu4fH9j/O+6X0KkgvYU72HuIg4vyLDhuISEB5U+M6BTqo7q7k159ag56qYGE40n2Be4rygtzOUlr4WbNI2TEBE66NZkbaCj6o+4gclPxj1DWoi09bfhnnA7DAtxs4mWh/t8Bl5OM3KjTxUPghwhMt6clKbB8xUdlSyKWvTqPcxLXIasRGxQWsQ5a3lPFv67KDsfneunXMtN2ffPJopjoiqQfjJoUuH0Ald0IkzgVCYXEh7f7vfvXNfOv0Svzr2K8wD5qD36U+bz1AQb4hnRdoKdpl2ucpgXD3ranSawJ9VYvQxGLQGjyYmJa4+GP+Dwpy4OcRFxHG8KbSRTEqI61ATEzjMTHXddSM6bSc7SgRTZnwmGqFhQeICr8esmIJGW4fJHaUvxFAU024oElWFEI6oxCA1iKePPs1ndZ9R11Xn8WXuD/737i+qBuEn++v3UzStaExqJBWmFAKOi3UkO6iU0lUY8EL7haBv8KVNpcRGxAbkUA+WjZkbeeSzR3jlzCu09rUGHN6qIIQg2ZjsUUAoN5tgciDct78wZWHIHdXuSXJDWTt7LRqhYXfV7lGZxyY6pg4T4OjCCA5B/kbFG9jstmE5OIozebQ+AXfSotNo6Wuh39Y/qHy4YilQfoOjJTshmw8ufuC3NUChoq2Cz2o/46HFD43KMjBaVA3CD9r72ilrKfO7k9pomZ84H71Gz6nmUyOOreqsct0gR2PrPNZ4jKIU/x3wo2Hd7HXoNDqeOfoMWqHlqln+Z08PxVuyXFlrGSmRKUFVoHVn4bSFVLRXBNxX2xcuAeFBg0gyJrFk+hI+qp7aTYRMZhMGrcEVmZSblEuvtZeqzqphY+u764nSRfmsXhAoij9jqKP6RPMJsuKzQravnIQcOvo7Aq46/NLplzBqjdyx4I6QzCNYVAHhBwcuHUAi/e7FPFr0Wj0LEhf45ahWGhMJRNC2zs6BTs63nw+47n2wxBviuTLtSnqsPSyZvmRUUVPJkck09w3XIMpay0ZlXlJYmLIQu7RzqmVkYe0vDd0N6DV6Eg2ey3ysz1jPubZzVJunbivSyo5KMuIyXNqCL0d1XZcjgimUPhlFMLmbmaSULgd1qAim5EZzbzNvX3ibm7Nv9loKZqxQBYQfHKg7QKw+NmRqpz8UpBRwuuW01yxhhSMNRxxtQJPzgtYgTjSfQCLD7qB2Z2PmRgC/mgP5IjlyeMG+fls/F9ov+N0a1RfujupQcalncJLcUK7JcJSSn8pahMlscpmXALLjs9FpdK7kRnfqu0OXA6GgaBDufr7arlpa+1pD6mcMprvcK+WvYLFbuDv/7pDNI1hUATECUkoO1B+gZEZJUI7UYClMKaTb0u2y1XrjcMNhlqUuG1XWZmljKQIxJg54hQ1zNvCthd/iluxbRrWdlMgU2vrasNqtrmUVbRXYpC0kGkSiMZHZsbNDmlHd0N3AjOgZXtfPipnFrJhZfpkYJyMWm4WazppB/i69Vs+8hHkeazLVdtWGLMRVYXrUdLRCO0iDUB4CQllJOdmYTLwh3u+Htz5rHzvO7GBt+lq/EnLDjSogRqC6s5rartox8z8oFCY7HdU+zEy1XbXUd9e7BERTb5NfzdiHUtpUSk5iTkCF8kaLUWfkO0u/Q4IxYVTbSTGmIJG09bW5lik3mVBoEODQIkJZk8lTktxQshOyqeiYvPHzvlCK9LlrEOAwM5W3lg+q1dU10EXnQGfIBYROoyM1KnVQuY3jTccxaA0hDWl2dZfz8+HtrfNv0d7fztaCrSGbw2hQBcQI7K/bDzBm/geFrPgsInWRPp8ilb7XxanFQWdtKhVqx9K8FEoUJ7R7JFN5azkx+hhmxYYm6qVoWhGNPY2u8NTRYJd2GnsaPUYwuZOdkI2pwzRIM5oqKDWYhkbM5Sbl0t7f7nLiQ3hyIBRmxswcpkHkJ+ej1wRXwdUbSqjrSC137dLOH0//kfzkfIpTi0M6h2C5bAWE0rBlJPbX72dm9EzmxM0J84wGo9VoyUvK86lBHL50mLiIOOYlzgs6a/NC+wU6LZ2TVkB4qsdU1lrGgqQFIYvICqUfoq2vDYvdMqKAyEnIwWK3eIzqmewMDXFVUDQ+dz9EOHIgFNKi01wCyGKzUNZSFlIHtUJ2QjadA5009fouHf9pzaeYzCa25m+dMEmSl6WA+Ov5v7LhtQ0jllCw2W0crD/IyrSV4/IPK0wp5EzrGSx2i8f1RxqOsDR1qaN0QPRMonRRAWsQSoLcZBcQShihzW7jXNu5kJmXwPFkq9foQ2Jm8hXi6s5Y1PHptfbylde+wocXPwzbPjxhMpuYFjltmElzfuJ8BGJQJFM4ciAU0mLSaOxpxGK3cLbtLAP2gbB0cvT34W376e2kRqWyIXNDyOcQLJelgLgy7Ur0Gj1/PP1Hn+NOtZyi09LJyplj639QKEwppN/W77E/RGNPI1WdVS5VNNj686VNpcQb4oc9zU0Whhbsu2i+SK+1N6RJZhHaCHKTckPiqHY1Cory7qQGmBs/F4EIax2fM61nuNR9iaONR8O2D09UdlR6dMBG6aOYEzdnkKO6vrueCE1EWPo6p8WkYZd2GrobXMJ/UUroH5T8EfanW05z6NIh7sq7K+QmrtEQVgEhhDAJIU4IIY4JIQ47lz0mhKh1LjsmhLjBbfyPhBAVQogzQoiN4ZpXcmQyN2XfxF/P/5XWvlav4xT/wxUzrwjXVHziy1Ht7n9QyE7IDviGcqzpGIumLZowKm2gROmjiNJFuUJdlSKAoYhgcmdhykJOtZwadWVXfzWISF0ks2JmhVWDUG7E1Z1jm28xNMTVHcVRrVDXVcfMmJlhSeBUQmfru+s50XSClMgUn9FlwZIcmUyiIdHn//Kl0y8RpYviq/O/GvL9j4ax0CDWSSkXSyndvS7/z7lssZTyHQAhRD5wJ1AAXAf8RggR2r6XbmzN30q/rZ8dZ3Z4HbO/fj95SXnjlqySHptOvCHeo6P68KXDROujWZC0wLUsJyGHlr4W2vva/dp+R38HlR2Vk9a8pOCeTV3eUk6EJiLkJUMWTltIr7V31E/0DT0N6ITOryficDecUW7ENZ01YdvHUNr62ujo7/D6/8lLyqO+u951Ddd314fFQQ1uuRBddSGr4OoNXw9vl7ov8X7l+3x13ldDmi0eCiaSiekW4BUpZb+UshKoAAIvheoncxPmsmrWKl4pf8VjH+AeSw+lTaWsSBvb6CV3hBAUJBd47A1xpOEIi6cvHpSbEWjWplKELpAS5hORlMgUVzZ1eWs58xLnhVxNV3JERuuovtR9ielR0/16Is5OyMZkNnn1QY0WxRlc21U7YoSNO5e6L42YwOkNV5E+LxqEovkp2k04ciAUFA2ivLUck9kUtkZZgMv8e6zx2LDXb0t/ix07d+XfFbb9B0u4BYQEdgkhjggh7nNb/pAQ4rgQ4gUhhPJ4Pgtw13VrnMvCxj0F99Da18rfLvxt2LrDDYex2q3j5n9QKEguoKK9gl5rr2tZa18r5zvODwuFCzRr81jTMTRCM6YZ4uEgOdJRsE9KGbISG0OZHTubBEPCqAWEPzkQCtkJ2VjtVqrMoY9kstgsVLRXpxZYWAAAIABJREFUEKuPpdfa63etoObeZq7feb3H34w/mMwmAK9JYEpwQXlrOX3WPlr7WkOeRa0QoY1gWuQ0Prj4ATD6DnK+yEvKo8vSxd3v3j3s9fq51/nKnK+ExRE/WsKdGnyVlLJOCDEd+EAIUQ78J/BTHMLjp8C/A/cCnnS7YY81TkFzH0BGRsawL1gsFmpqaujr6xtxcnHE8evCXyPNkrKywRmcEf0R/Ef+fxDfHk9Zx/DszrFilX4VBXkFnC0/S4Q2AnBkWz6d/zQpImXYvJ8peIbIvshhyz2xRC7hlwW/5GLFxbDM3RdGo5H09HT0+tE/6ScbkznQe4D67nrMA+aQRjApCCEoTCkcdenvhu4GvyvMuke/BNLX3B/Od5zHYrdwfdb1vHX+LWo6a/wqbFjRXoHVbuVo41Fuyr4p4P2aOkxEaCK8mo0SjAnMiJ5BWWuZK4ktXBoEOMJnjzcdR+DQ1sPFzTk3kxGXgcXmQRsUjGkVg0AIq4CQUtY5/zYKId4Alksp9yrrhRDPA287P9YA7k2J04FhBdullM8BzwEUFxcPEyA1NTXExsaSmZnplz1xZt9MartqSY8b3Jqyor0CndCNe7q7xeYIwZsRPcMV0lnfXY+hz0BuUu4wU4Whw4AGzYjzllJS3lpOvCE+rD9Ab/tuaWmhpqaGrKzR+wpSIlPoHOh03bxzk0OvQYDjR/xZ7Wd0W7qJ1kcH/H0pJQ09Daybvc6v8VnxWWiEJix+CMW8dG3Gtbx1/i2qO6v9KhV/scPxMFHeEly/iqFF+jyRm5RLeWu5KwcinNdnWnQax5uOk52QHdZKAnqNnpIZJWHbfrgIm4lJCBEthIhV3gMbgJNCCHd98TZAMbC/BdwphDAIIbJw9JY6GOh++/r6SE5O9tvZFGeIQ6fRDSr4ZrFZ6Lf2j2npCW/otXp0Gt0gE1OPpYcoXZRHO7ZBa6DPNrL21G/rxy7tROnC399iKEIIkpOT/dLy/EF58t1Xuw+N0DA/cX5ItjuUhdMWIpFB10jq6O+g39bvt4nJqDOSHpMellDX8tZyInWRjhwfhN+OasVEdK79XFBZ3iazacQAgrykPEwdJs53OARjuJzU8KXwCad5aTITTh9EKrBPCFGK40b/Nynle8BTztDX48A64H8BSClPAX8BTgPvAQ9K6aXX3ggEEomgERqSjEl0W7rpszpuWN2WboCgnhLDQaQu0iUgrHYrfdY+r42LDFoDNrttxB9vj9XR32AsGiB5IpTRIopm9VndZ2TFOUqUhAPlJhJswpwS4hpIKGUwuS3+UN5azoLEBRh1RlKjU6npCkxA9Nv6XQ5nf7HYHUX6Rsq5yUvKQyL5pPoTtELLtKhpAe0nEBThE44EualA2ASElPKClHKR81UgpXzCufxuKeVCKWWRlPJmKWW923eekFJmSykXSCnfDdfchpJoTEQjNC4tosvShVajxag1jtUUvLJ27VrOHD/DgG2A66+/nvpmx+nyJryUOXuKzHKn29KNTqObUEk5weJejylc5iVw9LHIT87n2dJn2XluZ0CRP+C7k5w3chJyqDJXebZdB4ld2ilvLXc589Nj0v3OhbhovujyjQTaFrW6sxqrtI6sQTiTHA83HCY1KjWsVZQLUgqI0ERwxYzxyXWa6EykMNdxQ6fRkWBIoGOgA4vN4rIxT5TkMaUl4qtvvoo+Wo8QwutTsjJW0YY8YbPb6BzoJC4ibsIc42hwd66Gw0Htzq/X/5pF0xfx6OeP8qN9P3Jpm/6g1P8KREBkJ2RjlVbXk3soqO6spsfa43KWp8em+2Vistgs1HbVsiZ9DQatwWNpbl94q8E0lNSoVBIMCdikLez+scKUQg5uOUhG3PCAFxVVQLhIjkxGSkl9dz1Wu5UYffD+B5PJRG5uLt/61rcoLCxky5YtfPjhh1x11VXMmzePgwcP0t3dzb333ktJSQlLlizhzTffBKC3t5c777yToqIivv71r9Pb2+uKXsqfl0/NpRoidZF89bavsmzZMgoKCnjuuedc+06MT+SZf32GVctXsWLFChoahlcg7bJ0IaUkzjCxknKCxT3pLBwhru6kRKbw7LXP8tDih3i38l3ufPtOv5+kG3oa0AptQG1Qg2k4MxLKjV05V7NjZ9PU2zTIz+WJ6q5q7NJOdkI28xPnB6xBjBTiqiCEcM1tLAIofDnML3fGrgPOOPB//nqK03Vmv8f32/qx2R0mnEh9B8JD5G1+WhyP3jRyOFxFRQWvvvoqzz33HCUlJfz5z39m3759vPXWW/zrv/4r+fn5XHPNNbzwwgu0t7ezfPlyrr32Wp599lmioqI4fvw4x48fZ+nSpQ5TkFaPHbvD/6CL4oUXXiApKYne3l5KSkq4/fbbSU5Opru7m+LlxTz86MP85xP/yfPPP8+Pf/zjQXMz95vRaXTj4qAOBxHaCOIi4jAPmMMuIMBxQ7l/0f0sTV3Kw3sfZsvftvCDkh9wx4I7fGpkDd0NpESmBHRDyozPRCM0IXVUl7eUo9PoXMInPSYdgNrOWnISc7x+z10DyEvK493Kd5FS+q2FVnZUkhKZMiha0Bt5SXkcqD8QthwIFf9QNQg3FHu8EBqPwiEQsrKyWLhwIRqNhoKCAtavX48QgoULF2Iymdi1axdPPvkkixcvZu3atfT19VFVVcXevXu56y5HRmVRURFFRY746EhdpMvmHa2P5plnnmHRokWsWLGC6upqzp07B0BERAQ3bLqBfls/S5cuxWQyDZqXzW6j0zJ1zEsKKZEppEWnjaq/daCUzCjh1ZtfpWRmCT/74mf8cO8PfdZqCiRJTsGgNZARmxFyDSInIQe91nG9p8c6BMRIjuqLZkeIa0ZcBrnJuXRaOv12boNDwPhbFHIsNQgV70xpDcKfJ313pJTUdtVi1BkDMgN4wmAwuN5rNBrXZ41Gg9VqRavV8vrrr7NgwYJh3/V043b3OXyx7ws+/PBD9u/fT1RUlEvAAOj1eow6I+397aABq3VwNNNUMy8prJu9znXDG0uSjEn8Zv1v+NXff8XzJ57ntnm3ee0+2NDT4HpqD4RgijB6Q8l/WZ2+2rVsdqwj/WgkR/VF80WSjEnEG+IHZTwr3x8Jk9nEtXOu9Wvs8pnLyU/OZ1nqMr/Gq4QHVYNwQwhBemz6qIWDP2zcuJFf/vKXLq3g73//OwCrV6/m5ZdfBuDkyZMcP+4IqVQEhFFnpKuzi8TERKKioigvL+fAgQODtq04qj2Funb0d0wp85LCd5d9lwcXPzgu+9YIDfcV3UekLpL3Te97HCOl5FL3pYAc1ArZCdlUd1YzYBsY7VRp7Gmkta91kCkuwZBAtD56REe1yWxyNc6alzgPrdAOau7ji7a+Ntr72/3WIFIiU9hx444xb9SlMhhVQIwTjzzyCBaLhaKiIgoLC3nkkUcA+Pa3v01XVxdFRUU89dRTLF/uqFeohK9G6aK47rrrsFqtFBUV8cgjj7BixeCCgoqAGFrkzWa30WXpIs4wtcxLEwGjzsja2WvZXbXbo2DusnTRa+0Nqpx0TkIONmkLOO/AE4pj2T3aSwjB7NjZfpmYlBu2QWsgKz7Lb0e14qAOdZVdlfAypU1M40VmZiYnT35ZgfXFF1/0uO7ZZ58d9t3IyEheeeUVj9s9f+E8eo0erUbLu+96ThPp6nKYkDRCw8abN/LNf/yma12npRMpJfERY2env5zYmLmRdyvf5eClg1yZduWgdUqjoGA1CHBEMrmXdw+GstYyBGLYdtJj0l2Zy57oGuiiubd5kAagOJL9QXFwZ8WpAmIyoWoQkwijzuhXBIwQAqPOOCxZToleClem8eXO1bOuJkoXxS7TrmHr/G0U5InMuEy0QuvTD2Gz23js88c4dOmQz22Vt5YzJ27OsETL9Nh0ajtrvZbxVhzU7gIiNymXpt6mQf3AvVFprkSv0atO50mGKiCmKAatgX5rv8vHoZiX4g3xqnkpTBi0BtZlrOPDqg+HmfeCyaJWiNBGkBHnO5Lpk5pPeP3c6/z74X/3meHtnkHtzuzY2QzYB2jsafT4PcVE5O4TUDKe/TEzVXZUMidujppzMMlQBcQUxaA1YJNf1mTqHHCYlyZax6qpxsY5G+no7+Bg/eA6kw3dDQhE0HWFchJyfJqAtp/ajkZoONVyymuP6Y7+Dmq7aj0KCCUXwpuj+qL5IgLB7LgvI5YUM5U/AiKQEFeViYMqIKYoiqNaMTOZB8zoNXrVvBRmrpx1JTH6mGHRTA09jiS5YGtfKZFMnmpsnWw+ydHGozy0+CESDAlsP7Xd4zY8OagVlFBVb45qk9lEWkya67oCiIuIY1bMrBEjmVxF+sa5dL5K4KgCYopi0H0pINTopbHDoDWwbvY6dlftHlRgL9gQV4XshGzs0u4xkmn7qe3E6GPYnLuZOxbcwZ7qPS6fgTuKgPBU0HBGzAw0QuM1F8KbBpCfnD+iBnG27axfRfpUJh6qgJii6IQOrUZLv61fNS+NMRszN2IeMA+K8Akmi9qdnPgvu8u5U9dVxwcXP+Af5v8DMREOIaHT6Pjj6T8O20ZZaxnTo6YPql2loNfomRk906OJSUo5KMTVndykXKo6q+ga6PI6951nd2LQGlg9a7XXMSoTE1VATFGEEK7mQR0DHSEzLz399NP09PS4Pt9www20t7cDEBMz/g2WJgIr01YSq48dZGZq6G4YlQYxJ24OOqEb5qh+ucyRVLklbwvgSDDbNHcTb1a8SXtf+6Cx5S3lPqvdpsemezQxNfc202Pt8SogwLsfomugi79e+CvXZV5HgjHBxxGqTERUATHJkVJit3sOTVQimbot3SExL9lstmEC4p133iEhQf3huxOhjWBdxjo+qvrIVT6+09I5Kg1Cr9UzJ27OIA2ic6CT18+9zobMDYMS8Lbmb6XP1serZ191Leu19lJprvRZzDA9xnPZb1cVVg8mJveSG554+8Lb9Fp7uTP3Tp/HpzIxUQVEiDGZTBQWFro+b9u2jcceewxwNP/57ne/y5VXXklhYSEHDzoiXR577DHuvvturrnmGubNm8fzzz/v+v4vfvELSkpKKCoq4tFHH3XtIy8vjwceeIClS5dSXT3Ybvzee++Rm5vLVzd8lZ89/DO+vfnbxEfE89hjj7Ft2zbXuMLCQlcxv1tvvdVj+fCYmBh+8pOfcMUVV/DEE09QV1fHunXrWLfO0Vc5MzOT5ubhcfCe5n05sTFzI52WTvbX7x9ViKs7Q7vL7Ty3k25LN/cU3DNo3LzEeVyVdhV/Lv+zqzzHubZz2KXdFZrqifTYdFr7Wof1uHDlQHhwMk+LmkayMdljbwgpJTvO7CA/OZ/ClMJh61UmPlM7k/rdh+HSidBuc8ZCuP7JoL/e3d3N559/zt69e7n33ntdWdXHjx/nwIEDdHd3s2TJEjZt2sTJkyc5d+4cBw8eRErJzTffzN69e8nIyODMmTP84Q9/4De/+c2g7ff19fFP//RPfPTRR8zMmMk/3PEPaIQGo853dzxf5cMLCwt5/PHHXeM+/vhjUlK816vatWuXx3mvXn352KBXzlxJbITDzHTj3BuB0QuInIQcPrj4Ab3WXvQaPS+XvUxxajEFycOLUm7N38r9H97PO5XvcGvOrT4jmBRckUydNYMyrU0dJiI0EV7LhOQm53rUII40HKGivYLHr3w8oONUmTioGsQYs3nzZsBRlM9sNrvs97fccguRkZGkpKSwbt06Dh48yK5du9i1axdLlixh6dKllJeXu8p6z5kzZ1gNJoDy8nKysrKYN28eRr2RG792IzqNbkTzkrfy4Vqtlttvvz2gY/Q178sFvVbPNbOv4aOqj1yRQaMxMYFDg5BIKjsq+eDiB9R317M1f6vHsSvTVpKTkMNLp19CSklZaxlxEXE++yu4yn4PMTNdNF8kIy4DjfB8u8hLyuNC+4VhxQR3nNlBbEQs12VdF8hhqkwgprYGMYon/WDR6XSDfAJKGW6FoTdq5bOn5VJKfvSjH3H//fcPWmcymYiO9tyT2n1bOo2O6VHTidBE+Jzbnj17vJYPNxqNaLWBZb96m/flxsbMjbx5/k3+u+K/AZgeNX1U23PvLvdy2ctkxmWyZvYaj2OFEGzN38pPPv8J++v3U9ZSRl5Sns8HBW+5ECazyWeZ8tykXKzSyrn2cy5tprm3mQ8vfsjmvM1q7s0kRtUgQkxqaiqNjY20tLTQ39/P22+/PWj9jh07ANi3bx/x8fHExzsK57355pv09fXR0tLCnj17KCkpYePGjbzwwgt0dTlCCGtra2ls9FwKQSE3N5fKykrOn3fYqt949Q3XTSEzM5OjRx1ZtkePHqWy0hFT39HR4bN8uDuxsbF0dnb6nEMw856KrJi5griIOE40nyDJmDQoySwYZsfNRqfR8drZ1zjVcoq78+/2+lQPsGnuJpKNybxw8gXOtZ0bsdteXEQccRFxg3IhrHbriEluLkd1y5dmpp3ndmKVVu6Yf4efR6cyEZnaGsQ4oNfrXU7drKwscnMH/ygTExO58sorMZvNvPDCC67ly5cvZ9OmTVRVVfHII4+QlpZGWloaZWVlrFzpaEATExPDn/70J59P9Eajkeeee45NmzaRkpLC1Vdf7fJz3H777bz00kssXryYkpIS5s+fD8B1113Hb3/7W4qKiliwYIFH05XCfffdx/XXX8/MmTP5+OOPPY7ZsGGDx3lPnz66J+jJhl6rZ33Get6oeGPU/gdw5CpkxmVytPEoCYYEbsq+yef4CG0Em3M386tjvwI8J8gNJT12cCRTXVcdVmn12ZchPTadGH2My1FttVt59eyrrJi5Qs2enuxIKSfta9myZXIop0+fHrZsorBmzRp56NChYcsfffRR+Ytf/CIs+/z444/lpk2bwrLt0TKR/1ehYl/NPln4YqF86MOHQrK9f9nzL7LwxUL5zNFn/Brf1tsmi/9YLAtfLJQVbRUjjv/enu/JTTu/vF4+qf5EFr5YKP/e8Hef37vn3Xvklr9tkVJKufvibln4YqH80PShX3NUGXuAw9KPe6xqYlJRCSPLZy5neuR05ibMDcn2ClMKidRFsjl3s1/jE4wJ3D7/dhINiX4Vy0uPSae2q9bVW1vp4zBSZ7e8pDzOtp3FZrex48wOpkdN9+ofUZk8qCamMWTPnj0elyt5EuFg7dq1rF27NmzbV/GNXqNn5y07Q+ao3ZK3hZuyb/JYLsMb3yv+HvcV3edXqe3ZsbOx2q009DSQFpPGRfNF4iLiSDD4TobMTcql19rLp7Wf8nnd5zy4+EF0GvX2MtlRNQgVlTATb4gnQhsRkm3pNLqAhAM4hJS/31FCXRVH9UXzRTLjM0cMk1Yc4E8efBKd0HH7vMBCo1UmJqqAUFFRcTE0F6LSXOmXaWpuwlwiNBHUdtWyfs76oPteqEwsVAGhoqLiYkbUDHRCR01XDT2WHhp7Gkf0P4BDS8lJdORKfH3B18M9TZUxQjUSqqiouNBqtKTFpFHdWU1VZxUwsoNaYU36GoxaI8WpxeGcosoYomoQYWAqlr1+8cUXqaurc33+1re+xenTpwHvBftUJidKLoSvKq6eeGDxA2y/frvalGoKoQqISYzNZhuz/QwVEL/73e/Iz88fk/2rjC2zY2dT01XDxQ5HFdeMuIxxnpHKeKEKiDAipeT73/8+hYWFLFy40FVm44EHHuCtt94C4LbbbuPee+8F4Pe//z0//vGPAfjTn/7E8uXLWbx4Mffff79LGLiX396/f/+g/R05coRFixaxcuVK137B8fT/0EMPucbdeOONrpDbb3/72xQXF1NQUDCoLHdmZiaPP/44V199Nf/1X//F4cOH2bJlC4sXL6a3t5e1a9dy+PDhYcfsbd4qk4f0mHQ6+js40XyCGdEz1FpKlzFT2gfx84M/H7FfbqDkJuXyw+U/9Gvszp07OXbsGKWlpTQ3N1NSUsLq1atZvXo1n376KTfffDO1tbXU19cDjvpMd955J2VlZezYsYPPPvsMvV7PAw88wMsvv8zWrVuHld9255vf/Ca//OUvWbNmDd///vf9muMTTzxBUlISNpuN9evXc/z4cYqKigBH2Y59+/YBDo1h27ZtFBd7ty/7mrfK5EEp2nfw0kGKphWN82xUxpMpLSDGm3379rF582a0Wi2pqamsWbOGQ4cOsWrVKp5++mlOnz5Nfn4+bW1t1NfXs3//fp555hm2b9/OkSNHKCkpAaC3t9dVx8hb+e2Ojg7a29tZs8aRvXr33Xfz7rvvjjjHv/zlLzz33HNYrVbq6+s5ffq0S0B8/euBRaPs3r3b67xVJg9KqGuvtddv/4PK1GRKCwh/n/TDhaPkyXBmzZpFW1sb7733HqtXr6a1tZW//OUvxMTEEBsbi5SSe+65h3/7t38b9l1v5bellF6dg97KfFdWVrJt2zYOHTpEYmIi3/jGNwaVJ/dVUtzb8Xqbt8rkQREQ4L+DWmVqovogwsjq1avZsWMHNpuNpqYm9u7dy/LlywFYuXIlTz/9NKtXr2bVqlVs27aNVatWAbB+/Xpee+01V4ns1tZWLl686HNfCQkJxMfHu0xCL7/8smtdZmYmx44dw263U11d7Wp1ajabiY6OJj4+noaGBp8ahz9lvoOZt8rEI1of7cq89jfEVWVqMqU1iPHmtttuY//+/SxatAghBE899RQzZjjaNq5atYpdu3aRk5PDnDlzaG1tdQmI/Px8fvazn7Fhwwbsdjt6vZ5f//rXzJnj+8f6hz/8gXvvvZeoqCg2btzoWn7VVVeRlZXFwoULKSwsZOnSpQAsWrSIJUuWUFBQwNy5c7nqqqu8bvsb3/gG//zP/0xkZOQw57hCsPNWmXikxzj6U6saxOWN8GYGmQwUFxfLoZE0ZWVl5OV577t7uWAymbjxxhtdvSAmIur/auLyw70/ZNfFXRzackgtujcFEUIckVKOmNGo/udVVFSGsTl3M4umLVKFw2VOWP/7QggT0AnYAKuUslgIkQTsADIBE3CHlLJNODys/wHcAPQA35BSHg3n/KYymZmZE1p7UJnYLJ6+mMXTF4/3NFTGmbFwUq+TUi52U2ceBnZLKecBu52fAa4H5jlf9wH/OQZzU1FRUVHxwnhEMd0CbHe+3w7c6rb8JWdHvANAghBiZjA7mMx+lcsF9X+kojLxCbeAkMAuIcQRIcR9zmWpUsp6AOdfJZNqFlDt9t0a57KAMBqNtLS0qDegCYyUkpaWFoxG43hPRUVFxQfh9kBdJaWsE0JMBz4QQviqe+Epy2vYXd4paO4DyMgYXkQsPT2dmpoampqagpyyylhgNBpJT08feaCKisq4EVYBIaWsc/5tFEK8ASwHGoQQM6WU9U4TUqNzeA0w2+3r6UAdQ5BSPgc8B44w16Hr9Xo9WVlZoT0QFRUVlcuQsJmYhBDRQohY5T2wATgJvAXc4xx2D/Cm8/1bwFbhYAXQoZiiVFRUVFTGnnBqEKnAG876QDrgz1LK94QQh4C/CCH+B1AFfM05/h0cIa4VOMJcvxnGuamoqKiojEDYBISU8gKwyMPyFmC9h+USeDBc81FRUVFRCYxJXWpDCNEEBFsNLgW43PtkqudAPQegnoPL8fjnSCmnjTRoUguI0SCEOOxPLZKpjHoO1HMA6jm43I/fF2q5bxUVFRUVj6gCQkVFRUXFI5ezgHhuvCcwAVDPgXoOQD0Hl/vxe+Wy9UGoqKioqPjmctYgVFRUVFR8MKUEhBDiBSFEoxDipNuyRUKI/UKIE0KIvwoh4tzWFTnXnXKuNzqXL3N+rhBCPOPsVTEpCOQcCCG2CCGOub3sQojFznWT8hwEePx6IcR25/IyIcSP3L5znRDijPP4H/a0r4lKgOcgQgjxB+fyUiHEWrfvTMprAEAIMVsI8bHz/3pKCPEd5/IkIcQHQohzzr+JzuXCeYwVQojjQoilbtu6xzn+nBDiHm/7nJJIKafMC1gNLAVOui07BKxxvr8X+KnzvQ44Dixyfk4GtM73B4GVOAoIvgtcP97HFo5zMOR7C4ELbp8n5TkI8Br4R+AV5/soHA2sMgEtcB6YC0QApUD+eB9bmM7Bg8AfnO+nA0cAzWS+Bpxznwksdb6PBc4C+cBTwMPO5Q8DP3e+v8F5jAJYAXzhXJ4EXHD+TXS+Txzv4xur15TSIKSUe4HWIYsXAHud7z8Abne+3wAcl1KWOr/bIqW0OQsIxkkp90vHFfISX/asmPAEeA7c2Qz8F8BkPgcBHr8EooUQOiASGADMOIpKVkgpL0gpB4BXcPQrmRQEeA7ycTTuQkrZCLQDxZP5GgBHKwHp7EgppewEynC0Dwi0H81G4AMpZauUsg3HubtuDA9lXJlSAsILJ4Gbne+/xpcVY+cDUgjxvhDiqBDiB87ls3BUllUIqi/FBMPbOXDn6zgFBFPvHHg7/teAbqAeR12wbVLKVkLUm2SC4e0clAK3CCF0QogsYJlz3ZS5BoQQmcAS4AsC70czFa8Fv7kcBMS9wINCiCM4VM0B53IdcDWwxfn3NiHEevzsSzHJ8HYOABBCXAH0SCkVm/VUOwfejn85jn7paUAW8D0hxFym3vGD93PwAo6b3mHgaeBzwMoUOQdCiBjgdeC7Ukqzr6Eelkkfyy8Lwt0waNyRUpbjMCchhJgPbHKuqgE+kVI2O9e9g8Nu+yccvSgUPPalmEz4OAcKd/Kl9gCOczNlzoGP4/9H4D0ppQVoFEJ8BhTjeGIcsTfJZMLbOZBSWoH/pYwTQnwOnAPamOTXgBBCj0M4vCyl3OlcHGg/mhpg7ZDle8I574nElNcghKObHUIIDfBj4LfOVe8DRUKIKKcNeg1w2ql2dgohVjijNrbyZc+KSYmPc6As+xoOOzvgUr2nzDnwcfxVwDXOCJZoHM7JchwO3XlCiCwhRAQOAfrW2M88dHg7B87rP9r5/iuAVUo56X8Hzjn/Hij7RmqGAAABBElEQVSTUv5ft1WB9qN5H9gghEh0RjxtcC67PBhvL3koXziegusBCw7J/z+A7+CIYDgLPIkzOdA5/i7gFA777FNuy4udy84Dv3L/zkR/BXEO1gIHPGxnUp6DQI4fiAFedV4Dp4Hvu23nBuf488D/Hu/jCuM5yATO4HDifoijyuekvgacc78ahynoOHDM+boBR7Tibhxa0m4gyTleAL92HusJoNhtW/fi6FNTAXxzvI9tLF9qJrWKioqKikemvIlJRUVFRSU4VAGhoqKiouIRVUCoqKioqHhEFRAqKioqKh5RBYSKioqKikdUAaGioqKi4hFVQKioqKioeEQVECoqKioqHvn/2roKHLgokj4AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x2aaad00e9470>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Plot\n",
"out=figdir+\"GB_transit_days_median.pdf\"\n",
"\n",
"years=np.arange(year1,year2+1)\n",
"plt.plot(years,median,label=\"median\")\n",
"plt.plot(years,uquartile,label=\"upper quartile\")\n",
"plt.plot(years,lquartile,label=\"lower quartile\")\n",
"plt.ylabel(\"Days\")\n",
"plt.legend()\n",
"plt.savefig(out,facecolor='w',format='pdf',bbox_inches='tight')"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"## Save to file\n",
"\n",
"dataout=xr.Dataset({'median': (['year'], median),'lquartile': (['year'], lquartile),'uquartile': (['year'],\n",
" uquartile)}, coords={'year': years})\n",
"dataout.to_netcdf(path=outfile)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SSH shkpwagn@hdata1.hlrn.de python_py3_std_/gfs2/work/shkpwagn/NB_WDIR",
"language": "",
"name": "rik_ssh_shkpwagn_hdata1_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
}
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"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\n",
"import itertools as it\n",
"from pandas import rolling_sum\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.6.4 | packaged by conda-forge | (default, Dec 23 2017, 16:31:06) \n",
"[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)]\n"
]
}
],
"source": [
"import sys\n",
"print(sys.version)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1948\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/gfs1/work/shkifmwr/_TM/software/miniconda3_20180131/envs/py3_std/bin/ipython:35: FutureWarning: pd.rolling_sum is deprecated for ndarrays and will be removed in a future version\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"1949\n",
"1950\n",
"1951\n",
"1952\n",
"1953\n",
"1954\n",
"1955\n",
"1956\n",
"1957\n",
"1958\n",
"1959\n",
"1960\n",
"1961\n",
"1962\n",
"1963\n",
"1964\n",
"1965\n",
"1966\n",
"1967\n",
"1968\n",
"1969\n",
"1970\n",
"1971\n",
"1972\n",
"1973\n",
"1974\n",
"1975\n",
"1976\n",
"1977\n",
"1978\n",
"1979\n",
"1980\n",
"1981\n",
"1982\n",
"1983\n",
"1984\n",
"1985\n",
"1986\n",
"1987\n",
"1988\n",
"1989\n",
"1990\n",
"1991\n",
"1992\n",
"1993\n",
"1994\n",
"1995\n",
"1996\n",
"1997\n",
"1998\n",
"1999\n",
"2000\n",
"2001\n",
"2002\n",
"2003\n",
"2004\n",
"2005\n",
"2006\n",
"2007\n"
]
}
],
"source": [
"## Load data\n",
"\n",
"wdir=\"/gfs2/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/\"\n",
"figdir=wdir+\"/FIGURES/\"\n",
"\n",
"\n",
"\n",
"fout15C=wdir+\"/DATA/15C_nfloats.nc\"\n",
"fout15C10d=wdir+\"/DATA/15C10d_nfloats.nc\"\n",
"fout10C=wdir+\"/DATA/10C_nfloats.nc\"\n",
"\n",
"\n",
"year1=1948\n",
"year2=1948\n",
"\n",
"n15C=[]\n",
"n15C10d=[]\n",
"n10C=[]\n",
"n15C10d_no10C\n",
"\n",
"\n",
"for year in range(year1,year2+1):\n",
" print(year)\n",
" ifile=wdir+\"/DATA/GS-\"+str(year)+\"/ariane_trajectories_qualitative.nc\"\n",
" ofile15C=wdir+\"/DATA/GS-\"+str(year)+\"/ariane_trajectories_qualitative_15C.nc\"\n",
" ofile10C=wdir+\"/DATA/GS-\"+str(year)+\"/ariane_trajectories_qualitative_10C.nc\"\n",
" ofile15C10d=wdir+\"/DATA/GS-\"+str(year)+\"/ariane_trajectories_qualitative_15C10d.nc\"\n",
"\n",
" data=xr.open_dataset(ifile)\n",
" #index15C=((data['traj_temp']<=15) & (data['nb_output'] < 73)).any('nb_output')\n",
" index10C=((data['traj_temp']<=10) & (data['nb_output'] < 73)).any('nb_output')\n",
" index15C10d=((rolling_sum(((data['traj_temp']<=15) & (data['nb_output'] < 73)).values,2)) == 2).any(axis=0) \n",
" index15C10d=xr.DataArray(index15C10d,dims=('ntraj'))\n",
" \n",
" index15C10d_not10C=(inedx15C10d & ~index10C)\n",
" \n",
" #data15C=data.where(index15C,drop=True)\n",
" #data10C=data.where(index10C,drop=True)\n",
" #data15C10d=data.where(index15C10d,drop=True)\n",
" #data15C10d_no10=data.where(index15C10d,drop=True)\n",
"\n",
" #n15C+=[index15C.sum().values]\n",
" #n15C10d+=[index15C10d.sum().values]\n",
" #n10C+=[index10C.sum().values]\n",
" #n15C10d_no10+=[index15C10d.sum().values]\n",
" \n",
" #Save to file\n",
" #data15C.to_netcdf(ofile15C)\n",
" #data15C10d.to_netcdf(ofile15C10d)\n",
" #data10C.to_netcdf(ofile10C)\n",
"\n",
"#ds15C = xr.Dataset({'ntraj': (['year'], n15C)},coords={'year':np.arange(year1,year2+1)})\n",
"#ds15C10d = xr.Dataset({'ntraj': (['year'], n15C10d)},coords={'year':np.arange(year1,year2+1)})\n",
"#ds10C = xr.Dataset({'ntraj': (['year'], n10C)},coords={'year':np.arange(year1,year2+1)})\n",
"\n",
"#ds15C.to_netcdf(fout15C)\n",
"#ds15C10d.to_netcdf(fout15C10d)\n",
"#ds10C.to_netcdf(fout10C)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "SSH shkpwagn@bdata2.hlrn.de python_py3_std_/gfs2/work/shkpwagn/NB_WDIR",
"language": "",
"name": "rik_ssh_shkpwagn_bdata2_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.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
{
"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\n",
"import itertools as it\n",
"from pandas import rolling_sum\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.6.4 | packaged by conda-forge | (default, Dec 23 2017, 16:31:06) \n",
"[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)]\n"
]
}
],
"source": [
"import sys\n",
"print(sys.version)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1948\n",
"1949\n",
"1950\n",
"1951\n",
"1952\n",
"1953\n",
"1954\n",
"1955\n",
"1956\n",
"1957\n",
"1958\n",
"1959\n",
"1960\n",
"1961\n",
"1962\n",
"1963\n",
"1964\n",
"1965\n",
"1966\n",
"1967\n",
"1968\n",
"1969\n",
"1970\n",
"1971\n",
"1972\n",
"1973\n",
"1974\n",
"1975\n",
"1976\n",
"1977\n",
"1978\n",
"1979\n",
"1980\n",
"1981\n",
"1982\n",
"1983\n",
"1984\n",
"1985\n",
"1986\n",
"1987\n",
"1988\n",
"1989\n",
"1990\n",
"1991\n",
"1992\n",
"1993\n",
"1994\n",
"1995\n",
"1996\n",
"1997\n",
"1998\n",
"1999\n",
"2000\n",
"2001\n",
"2002\n",
"2003\n",
"2004\n",
"2005\n",
"2006\n",
"2007\n"
]
}
],
"source": [
"## Load data\n",
"\n",
"wdir=\"/gfs2/work/shkpwagn/ARIANE/VIKING20-K301_Turtle/\"\n",
"fout=wdir+\"/DATA/10C_locations.nc\"\n",
"\n",
"year1=1948\n",
"year2=2007\n",
"\n",
"lon=[]\n",
"lat=[]\n",
"\n",
"for year in range(year1,year2+1):\n",
" print(year)\n",
" file=wdir+\"/DATA/GS-\"+str(year)+\"/ariane_trajectories_qualitative_10C.nc\"\n",
" data=xr.open_dataset(file)\n",
" ii=np.argmax(data['traj_temp']<=10,0)\n",
" lon=+data['traj_lon'].load()[ii,:].values\n",
" lat=+data['traj_lat'].load()[ii,:].values\n",
" \n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"lon1=int(lon.min())\n",
"lon2=int(lon.max()+1)\n",
"lat1=int(lat.min())\n",
"lat2=int(lat.max()+1)\n",
"\n",
"d=0.25\n",
"\n",
"xedges=np.arange(lon1,lon2+d,d)\n",
"yedges=np.arange(lat1,lat2+d,d)\n",
"x=xedges[:-1]+d\n",
"y=yedges[:-1]+d\n",
"\n",
"H,xedges,yedges=np.histogram2d(lon,lat,bins=(xedges, yedges))\n",
"ds = xr.Dataset({'hist': (['lat', 'lon'], H.T)},coords={'lat':y,'lon': x})\n",
"\n",
"ds.to_netcdf(fout)\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SSH shkpwagn@bdata2.hlrn.de python_py3_std_/gfs2/work/shkpwagn/NB_WDIR",
"language": "",
"name": "rik_ssh_shkpwagn_bdata2_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",