{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Tutorial on Units and UnitConverters in Parcels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In most applications, Parcels works with spherical meshes, where longitude and latitude are given in degrees, while depth is given in meters. But it is also possible to use flat meshes, where longitude and latitude are given in meters (note that the dimensions are then still called longitude and latitude for consistency reasons).\n", "\n", "In all cases, velocities are given in m/s. So Parcels seemlessly converts between meters and degrees, under the hood. For transparency, this tutorial explain how this works." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's first import the relevant modules, and create dictionaries for the U, V and temp data arrays, with the velocities 1 m/s and the temperature 20C. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from parcels import Field, FieldSet\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "xdim, ydim = (10, 20)\n", "data = {'U': np.ones((ydim, xdim), dtype=np.float32), \n", " 'V': np.ones((ydim, xdim), dtype=np.float32),\n", " 'temp': 20*np.ones((ydim, xdim), dtype=np.float32)}\n", "dims = {'lon': np.linspace(-15, 5, xdim, dtype=np.float32),\n", " 'lat': np.linspace(35, 60, ydim, dtype=np.float32)}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can convert these data and dims to a FieldSet object using FieldSet.from_data. We add the argument mesh='spherical' (this is the default option) to signal that all longitudes and latitudes are in degrees. \n", "\n", "Plotting the U field indeed shows a uniform 1m/s eastward flow. The .show() method recognises that this is a spherical mesh and hence plots the northwest European coastlines on top." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR4AAAETCAYAAADgV2kIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXeYHMW1t98zaWdzTsoSEggQCmSDyTkHm5wNxthgsE2wudfX9sc1ti/2BQe4ZBAy2WCSkU0WYLIASUgkZVhJm7V5d1Kf74/ulUarzbMzs7NT7/P0M9PV1VWnZ3fOVJ0KP1FVDAaDIZG4km2AwWBIP4zjMRgMCcc4HoPBkHCM4zEYDAnHOB6DwZBwjOMxGAwJxzieOCEiO4nIxyLSKiJXJtue0YqIHCwiVcm2w5BYjONxEBEVkek90n4lIg8Os8jrgEWqmquqfxaR+SLy69gt3cY+EZH/EZEG57hJRKSf/HNF5EMR6XBe546Gsgzph3E88WMysGKkChMRTy/JlwInA3OA2cDxwPf6uN8HPAM8CBQCDwDPOOlJKytROM7Q/L+PFlTVHPbsbQWm90j7FfBgH/l3AF4FGoB64CGgwLn2KhABuoA27C9iCAg65885+cYBTwJ1wFrgyh51P4H95W4BLunFhreBS6POLwbe7cPeI4ENgESlfQUcneSyDgaqgKuBWmATcFHU9XxggfMZrQd+Drh6+/sAU5y/o8c5XwTcCLwFdALTgQuBNUCr85mfk+z/vXQ8zC/A8BHgt9jOY2dgIvYXAVU9FHgTuEJVc1T1LmzHdJNzfoLz6/scsBQYDxwG/EhEjoqq4yRs51Pg3N+TXZ37u1nqpNkGiiwTkbOj8i5T5xvpsCwqf8LK6oUKbAczHttJ3SYihc61vzjXpgEHAecDF/VTVk/Ow3b8udjO68/AMaqaC+wHLBlCWYYRorfmu2EQqOoqYJVzWiciNwO/HEIRewGlqnqDc75GRO4GzgRecNLeUdWnnfedvZSRAzRHnTcDOSIiajO7n7zd+XMTXVYvzxECblDVMLBQRNqAnUTkA+AMYJ6qtgKtIvK/2M7k3l7K6Y35qroCQETCgAXMEpGvVHUTdgvLkGBMi2crEcDbI82L/aXYDhEpE5FHRWSDiLRgd4lKhlDfZGCciDR1H8B/AOVReb4eoIw2IC/qPA9o6+PL3TNvd/7WJJcF0OA4nW46sJ1XCeDD7mJ1sx67ZTRYtnyGqtqO7cguAzaJyPMiMnMIZRlGCON4tvIVdowgmqls+08fzW+x4wmzVTUPOBe7+9UXPb90XwNrVbUg6shV1WP7uacnK7ADuN3Moe+A9gpgdo/RpdlR+ZNVVn/UYzv+yVFpk7DjSwDtQFbUtYpeytjmM1TVF1T1CKAS+By4exh2GWLEOJ6tPAb8XEQmiIhLRA4HTsCOsfRGLvYve5OIjAeuHaD8Guw4RTfvAy0i8lMRyRQRt4jMEpG9hmDzAuAnIjJeRMZhB2jn95F3EXar7koRyRCRK5z0V5NcVp+oagR4HLhRRHJFZDLwE+zWJdjxmQNFZJKI5APX91eeiJSLyIkikg0EsP9+kaHaZRgBkh3dHi0HkAn8HliHHZP4CDixn/y7Ah9i//Muwf5yVUVdX0TUSBQww8nXBDztpI0DHgGqgc3Au8DhzrVf0ceIWlSZAtwENDrHTWw70rSCqFEbYJ5jc6fzfPOSUVaPZzg4+nNz0tZFfQ6F2I6mDruV+AucUS3n+m3OZ7oK+C7bj2pF/w0qgdedv2+Tc32XZP/vpeMhzh/EYDAYEobpahkMhoRjHI/BkEaIyH0iUisiy/u4PlNE3hGRgIhc0+Pa0SLyhYisEpGfRaVPFZH3RGSliDwWNYO9T4zjMRjSi/nA0f1cbwSuBP4QnSgibux42jHALsBZIrKLc/l/gFtUdQZ2rPLigYwwjsdgSCNU9Q1s59LX9VpV/YDt56/tDaxS1TWqGgQeBU5yplQcytbR3wew1+n1y6iduVxQWqrN9fXJNsNgiBfrVXVK98lRh2RrQ2NsI/sfLguswF4f2M1dai/XGQnGs+2E1ipgH6AYaNKtE0CrGMQEz1HreJrr65n6pz8MnHEE+EnFBG6uTtyWMGO5PvNsg2PtVddET4qkoTHC+y9MiqlMd+XKLlXdM6ZC+qa3ybHaT3q/jFrHYzCkEwpYWMk2oz+qsBdCdzMB2Ig9u7xARDxOq6c7vV9MjMdgGBUoEbViOuLMB8AMZwTLh72Y+Vm1JwK+BnzbyXcB9l5N/WJaPAbDKECBcAJWb4jII9izxUvE3nL2lziLo1X1DhGpABZjL+y1RORH2LO7W5ylMS8AbuA+dVb9Az8FHhV7h82PGcTOAcbxGAyjAEWJJGAVgaqeNcD1auzuUm/XFgILe0lfgz3qNWiM4zEYRgnWwDHZMYNxPAbDKECBiHE8BoMh0ZgWj8FgSCgKCYnxjBaM4zEYRgmjehbPCGMcj8EwClDUxHgM6YOq0v7hx1jBIJ68PNz5ebjz8nDn5iCu7eeXqmVhdQWwOjqco5PIlvcdWGee3UsthgFRiKSP3zGOJ91p+ucLtC/9hIzJk+loaSHS0kK4uQWrowN3Tg7uvFwQwerqsh1LZxfi8+HOysSVlYUrKwt3dhbhlhYCa9fD6Wcm+5FSEkUI9asVMLYwjieNaXr5VdqXLKPyh9/HnZu7zTWNRIi0thJpbUMjEdxZWbiyMnFlZiJu9zZ5gxs3UX37XZSedzYuvz+RjzBmUMAyLR7DWKf5jX/T+s57VP7wB9s5HQBxu/EUFOApKOi3nGB1NdW3303RySeSs/vceJmbFkRMi8cwlml95z2aX11E5Q9/gKcgf9jlBGtqqf6/uyg68Thy9pg3ghamH/YEQuN4DGOUtsUfsfmfL1B5xffxFhcNu5xgTS3Vt91J4fHHkrPXHiNoYfpiqXE8hjFI+9JPaHzmOSq+fynestJhlxOqraP6/+6k8NijyN07XvtOpRemxWMYk3R89jn1jz9BxWXfxTeuctjlhOrq2XTbnRQefSS5+w5pQbKhHxQhkkbbYxnHkwaoZVG34GHKL/0OGRN73fFgUITqG9h02x0UHHU4ud/YZwQtNIDpahnGGFZHJ5k7zcA/dcqwywg1NFJ92x0UHHYIufvsRaihkXB9A6GGBsL1DYQbGgg1biZ49TVs+P0t9k1if5HEJYDY5yKI20XJOWfiLRp+jGmsYbpahjFHpKOD7BhGnayuLqpvvZ1IWzvNi96k8enncOXk4C0pwlNcgrekiKzZu+EpLsJTXETJmac5N9qrj1QVnAWQVleA2nvux+XPjPm5xhKKENL0+Tqmz5OmKeHNTWheEVk7zxx+IS4XBUcejjsvD09JMZ6iQlxeb+9Zvd5+u3PtHy/FP30H3FnG8fTEtHgMY4b2j5fgOvo4JDj8tc8un2/EYjrtyz4ha85uI1LWWEJViGj6BJcH9aQisk5EPhGRJSKy2El7zDlf4lxfEpX/9yKyWEQOcs6niIiKyA+j8twqIheO8PMYetD24ce4s7OSbQYAVihE52efkz1r12SbMiqxkJiOVGIoLvYQVZ3bLRimqmc453OBJ4G/gy367uQ/ELg86v5a4KrBCLobRoZgdQ2R1lYkIyPZpgDQ9eVKvOMq7YWnhm2wg8uumI7BICL3iUitiCzv47qIyJ9FZJWILBOR3Z30Q6IaGktEpEtETnauzReRtVHXBlw7E3NXy9FOPh1bPxls6QuL7VUG64C3sHV37o61XsPAtH/0Mdnz5o6a38L2pcvJnm26Wb2TsK7WfOBWYEEf148BZjjHPsDtwD6q+howF0BEioBVwItR912rqk8wSAbreBR4UUQUuLOHHvMBQI2qrgRQ1RUikgX8G7i2Rzm/A/4pIvcNptKfVAx/zslQKPf6ElZXoupTIHTM8XiKixP6fH3VpUDwjDPxVZRvt7o9HvXFg5Gs64c9zm0l0fg7HlV9Q0Sm9JPlJGCBI9T3rogUiEilqm6KyvNt4J+q2jFcOwbrePZX1Y0iUga8JCKfq+obzrWzgEeiM6tqz8+1O32tiLwPDGq3qFTUxB4t9QXWf0XtgoeZ8POfcrVvYlI/S7Us6h9/kkhTExWXfTfu9cWLeNcViX0CYUl3DNbhrh6NhMEwHvg66rzKSYt2PGcCN/e470YR+QXwCvAzVQ30V8mgHI+qbnRea0XkKWzxrjdExAOcCgxlleBvgCeANwbKaBg+rpwcNBCgY+knUDlx4BviRLfTCdXWUvG9S5Jmx2hnhJZM1HfHYGOgN++3ZacgEakEdsNWFO3meqAa8AF3YSuL3tBfJQM+qYhki0hu93vgSKA7MHU48LmqDvpnQFU/Bz4Fjh/sPYah4y0uovyyS2h44ikinZ1JsUEti/rHntjidFyjJMg9GlEgpJ6YjhGiCoj+pZoAbIw6Px14SlVDW2xX3aQ2AeB+BqEqOhgXWw78W0SWAu8Dz6vqv5xrZ9KjmzVIbqQPmVTDyJExYTzl37uYcONm2j/pdRAjbqhlUf/oE4Tr6o3TGQSKENHYjhHiWeB8Z3RrX6C5R3xnu9CK0wrqHmg6ma0Nkz4Z0E06ushz+rh24UD3O/nWAbOizpcytKF8wzDJmDgBb1Ep9Y8+gYiQlYA5NFucTn095d+72DidQZKI4LKIPAIcjB0PqgJ+CXgBVPUObG30Y7FHrTqAi6LunYLdGnq9R7EPiUgpdjdtCXDZQHaYmctpgMvnpeLS71B9171UXn5ZTNtiDIQC9Y/+jXBDg3E6Q0CVhAynq+pZA1xXtp1/F31tHXaguWf6odvn7h/T6kgTMiZPInP6DgRrauJWh1oW4cZGwg2NlF9qnM7QiG3W8lieuWxIcSQzE6uzKy5lq2VR/8jjaDhC+aXfMU5niNgSxq6YjlTCdLXSCHemHytOI1ydX66k7f3F6Jlns/EPf8SV6Uf8mbgz/bgyM7cc7rxccvbaY0QnEY4VzA6EhjGJ+Ldv8WgkQripCW9xcUxlZ83cicm//y2+inLKLrkIq7MTq7PLee3cct7+1nJCtXUUnXhcTPWNNRQxOxAaxibuTD/BlpYt51YgQO38v9L15SqKv31KzFtfuHxexO3GV17WZ55IWzsbfn8L/mlTyZq1S0z1jTXSqcWTPk9q2CbGE25pYdNfbsedm8e4a35E0yuv0fDUs6g1/H17BoM7J5uyC86h7pHHCTdujmtdqYS9A6E7piOVMI4njeiO8QRratn0x1vJmrULJWedhq+ygnE/vpLA119T/9igFxgPG/+0qeQfehC1DzyIhsNxry8VsCWMXTEdqURqWWuICfH7CW7cyKY//x8FRx1B4dFHIlEbslutbfinTE6ILfmHHIQrO4vG5xYmpL5UIILEdKQSxvGkEZ68PKyuAKXnnUXuPnttSVfLonbBw/h3mpEw2RpxuSg95yw6ln2S8OUcoxFVSasWjwkupwCR9g6CVRuwAgE0EMByDvt9EA0GyZ47h8ydZvRbjreslMm/+X/bDWVvfv5faDBA8SknxfMxtsOdnUXpBedSc/d9+MZVxjyyluqk2lycWDCOZxSjlkXb+4tpfG4h3vJS3JmZiC8Dlz8DycjAleHDk58HItQ/9jd848ZRdPLxeEtK+iyzp9PpWL6Ctg8/Yvw1P0rK3Br/lMkUHHEYtff/lXE/ugLxpOe/pL0RWGp1l2IhPf/KKUBwUzX1f3sSDYaouOySARVAc7+5Hy2L3mDjzX8md999KDjyMFx+/4D1uPPy0GCIcONm3Dk5I2X+kMg76AC6Vq+h8Zl/UPytk5NiQ/IxKhOGJGIFgzQ+t5BNf7mdnHlzGfeTKwclO+zyeik44jDG//RqIq2tVN14E63vfTDg8HjGpImUnHUaNXffT6i2bqQeY0iICCVnnUHHik9pX7I0KTYkG3tUS2I6UgnT4hlFdKz4jIYn/k7GlMmM/+nVdjdqiHjy8yk950wC67+i4e/P0PLmW1i/+U2/92TvNotISysbfn8z4svAW1riHKV4y0rwlJbiLSmO6/ord1YmZReeR/Wd9+CbML7f7uJYJZ0mEBrHMwoI1dXT+MxzBDdVU3zGt8mauVPMZWZMnkTlj66g/cOPCdU3UPvXhyg68Tg8hQW95s/b/xvk7rcvkeYWQnV1hOrqCdfV07b4Y0J1dYQbGnBlZeEttR2Rf/Ikcvbde8tw/EiQMWkihUcdQc09D1Bw+CH4p03FU1Q4YuWPZhQhnGKTAGPBOJ4kYnV10fTiy7S++z75hxxE6QXn9ikNPBxEhJw9d8dXXoE7J5uq39zEpF//ss+Wi4jgKcjHU5BP5ozp21xTyyLS1Eyorp5QXR0tb71DYP1XFJ/+LcQ1cr/UuQfsDx4P7cs+oeGpZxGPB/+0qfinTcE/bSrecZUj6uxGC/Z+PGPvufrCOJ4koEDru++z+fl/kTlzR8b/9JphdasGjQihunpy99t32N0lcbnwFBXiKSokc6cZ5Oy5OzX3PkDt/X+l9PyzR8xhigh5++1L3n77oqqE6xvoWrOWrtVraXr5VfIPOoD8Qw8ekbpGG6kWp4mFYUsYR127xpEnLnHOXSKyQETeFpFdnbSDnTwnRN33DxE5eASfJSXoWr2WUE0Nre++T/klF1J6zpnxdTpApLkFDQZHdEW4y++n4nsXg0uoufNerK6R3+dHRPCWlpC7z16Unn06FT/4Hk0vv0akfdhyTqMWe3V6+kwgHLaEMYCITASOAL6Kynck8B5wCnB1VHoV8J8x2JrShBs3Uzv/QWoXPIg7N5fKqy4nY/KkuNfb8elnRNrbKb3w3BGfpyMeD2UXnIu3tIRNt95BpK1tRMvvia+8jOw5u9H88qtxrSdZmCUTg+cW4DqidHfYKmFssa1Gz1KgWUSOiLHOlMIKBNi88AU2/OEWvOWlTPiP63BnZSUkThFqaKDuocfwFhfjyYtPq0pcLopP/xZZO89k059uQyORuNTTTcHRR9D67vtjbmV7ug2nD9bxdEsYfygilwKIyInABkcxIpoXgIOwZTJ6qg3+Gvh5DPamDKpK2+KPqPrNTYTq6hh/zY8pPOaohG0JagVD1N63gIIjDsOV4YtrXSJC4XFHk7v/fgRrauPS7erGk59P7je/weZ/vjBw5pQiMV0tEblPRGpFpNcFco6szZ9FZJWILBOR3aOuRZxwyxIReTYqfaqIvCciK0XkMREZ8B9u2BLG2N2mI3tmVNUwtt7WdqjqmyKCiBwwmEpTURO7m1BdPbrXPniOOHq7L34i9L7DzS3o9f+BpyRx2ul6xpmUIVyNxrVFpxddTHBTNd7i0jGjnQ4JWzIxH7gVWNDH9WOAGc6xD3C78wrQqapze7nnf4BbVPVREbkDuNi5r0+GK2F8EDAVWOr8g00APhKRvVW1eoDibsR2WgNuxJLKmtg19z6Af+rkXkdgEqH3XX3XvWTP3o3cffdOmL545+dfcu0us7jVCsa9ruY336Dzi5Xc9MtfpvT/STeJGk5X1Tccfay+OAlY4MjcvCsiBSJS2UPUbwuOiN+hwNlO0gPArxjA8QxXwvgDVS1T1SmqOgU7cLz7IJwOqvoiUEgfIoFjhcJjjqTplUVx7Xb0h3+HaQSrNiS0zs6VqxLWlcz75n6EqquxAoGE1BdvFCFsuWM6sEX6Fkcdlw7DlPHA11HnVWzV0vI75b4rIt2L6oqBJqen0zN/nwymxVMOPOW0bDzAw1ESxsPlRuCZGMsY1fjGVZK54wxa3vg3BUcenvD6/dOm0vDEU8O6N9LailpqK0V4vYPuNnWtWo0rww+drcOqdyiIx0PhsUcTbmpGNb5du0QxAl2t+uhR52HSmxHdg0eTnJDLNOBVEfkEaOknf5/EJGEclWfKANcXAYuizp+l9wccUxQecyQb/3grud/cH3dWZkLrzpg4gVDt0AO9GolQ9ZvfI243VmcnqorL3y1RE/Xq7z530jIyCG7ciMQ5kB1N9h7zQJWOZcvJnrNbwuqNB92jWqOAKmyZ4m4mAN2hlu7XNSKyCJgHPAkUiIjHafVsyd8fZuZyHPGWlZI1axdaFr1O4bFHJ7Ru8XjwTZhAYP1XMGX6wDc4dK1Zi6ekmPFXXwWAFQqhnV1YXV3bStZEnYfr6rG6usg/9OCEtjzE5cJTkM/mfywka9edU34vn1EyCfBZ4AoReRQ7qNysqptEpBDoUNWAM1l4f+AmVVUReQ34NvAocAGD6M2k9l8qBSg46nA2/uFP5B14AO6c7ITW7Z86ha41a+2hgEHS+elnZO0yc8u5y+sFrxd3Xu7IGzgCiN+Pt6Kcxmefp/jUxO6gOBxUlVB1LzLSCZqLIyKPAAdjx4OqgF8CXse2O4CFwLHAKqADuMi5dWfgThGxsGPDv1PVT51rPwUeFZFfAx8D9w5kh3E8ccZbXEz2vDk0v/oaRScen9jK3S7CDY1DuqXj088pPfuMOBk08ghQcuZpbLjpZjJn7kjWLjsn26ReibS20vTya3QsX9HrJMtE7UCoqmcNcF2By3tJfxvotT/rhGP2Hoodo6JtN9YpOPIwWt95j3BLb3G4+NCxfAVt735AwbFHDfqeUEMjkbY2fIPYeGw04c7OpvTcs6h/5G9EWuIf2B4ODU8+TaS1jbLvXMDEX/a+ciidZi6bFk8C8BQUkLPXnjS//FpCugOBqg3UPfw45Zd+B29R0aDv6/zsc7J2njmi21wkiswZ08n9xj5U/e4P9sbx5WV4y8vwOa/u/PykjXx1frmSwPqvGH/9tbh8vQffR1FwOSEYx5MoVLE6O+NeTbi5hZq776f4tFOGrJEVbmom0t6ORiJJ2fg9VgqOOZLcffcmWFtLqLqWUHUNHUs/sUf3AkG8ZaVbnVFZGRlTJuEp6H1jtJFCIxEannyaopNP7NPpdGMcj2FEaft4CR0rPmPcNVfFtR4rGKTm7vvI3W9fcub1NrO9fwqPPsLeY2f+Xym74NyUGyUSkS17BtFjF8dIRweh2jpCNbWEamppXvQG4nJReeUP4mpTyxv/xlOQT9bsWf3ms3cgTL2W5nBJnydNEsHqahr+9hRl3zkfd1ZW3OpRy6LuwUfwlpdRcORhwypDPB7KL74AVKm9f8GYkhd2Z2XhnzKZ3H32oujE48jefS7e8rK41hlubqHppVcoOvXkgbt5ml4xHuN44ojV1UXtfQsoOvE4MiYMOIs8JjY//y8iLa2UnnV6TLEM8Xgou/A8cLmpuXc+Vig0glaOHoJVG/CNHxfXOjY/9zy5++6DbxAOzmyLYRgRVJW6Rx7HP20qufsOaaRxyLS+/wHtHy+l/JILR6R7ZDufc3H5fNTeMx8rOPacT7BqA744/hh0rV5L58pVFBw1+OUy6eR4UqsTn0K0LHqTcEMjpVdtNyViRFBVAuu/ov3jpbQt/ojKH35/RAT5VJVQTQ0dKz4j3NxCcMMGglVV+KdNHQGrRwcaDhOqrcU3rjI+5UciNDzxd4pOOmHQi2btrU9Ty3nEgnE8caBr9RqaX3mNyp/8cERVI6KdTfvSZbh8PrLnzmbcVZfjLSsdfrnhMJ2rVtO54jM6VnyKWhZZu+5CwZGH4Z8+HZdv5J5hNBDcVI2nqGjAUabh0vr2u7iys8ieN7QNGNQ4HsNwCTe3UPvAg5Scc8aQ5tD0RV/OpuLSi/FWVgw7nqORCG0ffEjH8k/pXLkKX2UFWbvuTPklF8VU7mhGLYvWd96j6Z8vknfIgXGpI9LWxuZ/vUjlFd8f8mdotNMNw0IjEWof+Cu539iXrJ1nDnxDX+XEydlEE2lrp+Fvfyf/sEMoOfO0hK8jSzSdK1fR+NSzuPx+yi+7JG7B/sbnFpKz5x74KiuGdJ+qmcdjGCbNi94gsO4rMiZMoPm113EX5OPJz7df8/L6Dfz27mzmjKizicaTn0fhScfT9sGHQwqApiK1Cx4msHYtRScdT9ac2XFrzQXWf0Xnp58z4T+vG9b9pqtlGBY5e8zDk59PuLmZcONmutauI9LUTLi5mUhrG67MTDwF+YSuv576x5/EU5CPOy+P4KbqhDibnuQdsD+dn35O079eovC4xG7bkUg0EiZz55lkz43fppdqWdT/7e8UnngcLr9/GCUIESt9BpmN4xlBPAUF5Oy5e6/X1LKItLYSaWrGlZ2Nr7KCcFMzoVWr8RQVJczZRCMilJ59Bht+fzORzg48hYVbW2gF+bjz80c0OJ4sSs88jQ3/+yfaPviQnL32GPHyw5ubaHxuIeL19vn3HwizVssQF8TlwpPvdL0yM8k7YP9kmwSAOy+Xiu9fSucXKwk3NRH8eoPdQmtqItzcgisjY4sT6u462ud5eAoK8BQXJWyf5eHiysyk/KLz2XTbnfjGjxuxYfRIeztNL71K23sfkLv/Nyg57ZTh/3CoHedJF4zjMeAbV9nrl1EtC6ujw1482tRMuKmJzi++pPnlV7csp8iasxvl37kg0SYPGd/4cRSdciK19y9g3NVXDbM7ZMfiwnX1tH28hJbX3yR77hzG/+yaEZGhNqNaPRCRdUArEAHCqrqniJyGLWOxM7C3qi6Oyv974BDgalV93ZHTWAtcqap/cfLcCixW1fkj9TCGkSfc0EjX6jV0rVlL1+q1uLKyyN5jd/w7TMW/wzQ8xbFPGUgU2bNn0fzyqzS98BJFJ50wqHu6J1R2rVpD6JDD+foX/w0uIXPmToz78ZV4S0tGxDbFBJf74hBVrY86Xw6cCtwZnUlEuseRD8QWD3vdOa8FrhKRO1U1/sJLhpgINTRS99eHsTo68M+YTvbu8yg+7VQ8+fnJNm3IqCrtHy2h8dl/4J82jbyD+p7Do5EIwY2b6Fq9lq41a+havQaXz4d/+g64/BlUXnU5nuKiOMTizMzlQaGqnwG9/QG6tdOVbZUk6oC3sDeDvnu49RriT9vij2h46hkKDj+UvIMOSMmNwboJVG2g4cmn0WCQsvPPxb/Dtks/rM5OutZ9RWDtWrrWrCPw1dd4CgvwT5tK1m6zKD75RHubDeydDr0lxXGz1cR4tqdbO12BO1X1rj4zqq4QkSzg38C1PS7/DviniNw3mEpTUZo2letTVcKNm9E998Zz1DHDGtEaLc+mlkW4qRmrciKe//w5rpzsbX71omOWAAAgAElEQVQFI52dRJpbUK8gu87GtfueuDIyEJ+vT0cbbwlj09Xanu2001X1jb4yq2pvnyuqulZE3mer3Gm/jAVp2lSpr2vdeuoWPEzmTjMoOuVEXA3D24tnpJ8t0tEJlgUiIIK4nC+nCIiLqydO4eaNX21NU6XlzbdoevEVcvbcnYKjj8TdZkHb5i1lqipVv/4dRSccR9asXRBLoTMwoBBhvCWME+F4nB/944FaVd1udzJHkvhP2EoTHcCFqvqRiMzFliXOw4713qiqjzn3zMfWMml2irlQVZf0Z8dwtdP3Bvp0PAPwG+CJGO43jCBqWTS99Aqtb75N8enfInuAnfL6LUsVKxAg0tE5JAFDVSXS1EywpoZQdQ2hmtot7zUcQTzu7m8mamm34aBK8Le/Y91113cXBKpkztyJyit/gK+ivNf6AuvWI243WXN2G1Vr0iJWQmyZD9wKLOjj+jHADOfYB9vZ7IPthM5X1ZUiMg74UEReUNUm575rVfWJwRoxoONx9NJdqtoapZ1+w2Ar6Imqfi4in2J73feHW44hdsKNm6l98GHE5WbcNT/CUzD8wHG4pYX6hx8j/L3vU3XDb8ieN4e8A/bfZpheIxHCDY2OU6klVFND0NmK1JXhw1tejreiDN+4yi07BLpzc/t1DhkVE5h6y01b6xiEnHHb+4vJ2WuPUeV0IDEtHlV9wxll7ouTgAWOzM27IlIgIpWq+mVUGRtFpBYoBZr6Kqg/hq2dLiKnAH9xKn9eRJao6mC1VG7EFv4yJIm2j5bQ8ORT5B96MPmHHBRTALn9kxU0PP4Eud/YB29lBeOvv5bWt9+l+va78ZaW4M7LJVhdQ7i+HnduLt6Kcrzl5fh3mEbu/t+wHcwIbQs7kDPRcJj2JcsYf+2PR6S+kUKRkXA8JSKyOOr8rv7isX0wHvg66rzKSdvUnSAiewM+YHVUvhtF5BfAK8DPVDXQXyXD1k5X1aeApwa638m7DpgVdb4Us/thUrC6umh48mkC69ZTcdl3yehFQ0sjEdSyBgwuW4EAjU8/R+fnX1B20fn4p01FcBagHnMkBUccSsfyT9FwmPzDD8VbVhq3PXAGS8fyT/GNr9wyUjWaGIFBrXpV3TPGMnrzfltME5FK4K/ABapqOcnXA9XYzugubGXRfntFZuZyGhFY/xW1Cx7CP2M646750XZLHdSyaHt/MZv/+QKR1jbcWVl4igrwFBXZ67iKuo8irK4u6h9+jIzJkxh/3U9wZW4f0xGPh+y5sxP1eIPCXq8V63czDiQouDwIqoCJUecTgI0AIpIHPA/8XFXf7c6gqt2toYCI3A9cM1AlxvGkCc2L3qTpxZcpOf1b2zkDVaXz089pfO55XH4/ZReeR8bkSURaWwk3bnaORoIbN9Gx/FPCmzdjBQIUnXQ8ObvPS9ITDZ1IWxtdq1dTel6/Kr7JY3TM43kWuEJEHsUOKjer6iYR8WH3cBao6t+ib3BiQJucEbGTsScX94txPGmC1dFuDzW//S4aDpM1ezdcPi+B9V/R+OzzRFpbKTzhWLJm7bolTtK9qJWpU5JqeyyoZdG1eo2zz9En5Oy917DXacWbBA2nPwIcjB0PqgJ+CXjt+vUOYCH2UPoq7JGsi5xbT8dejVAsIhc6ad3D5g+JSCl2N20JcNlAdhjHkyYUHns0+YcfRsfy5bS9t5iGJ5/CN24cobo6Co4+ktx99kpJ9dDe6OlsPAX5ZM+dw7gfX4G3ZGTWVsWDRMxcVtV+m3vOaNZ2CgWq+iDwYB/3HDpUO4zjSSNcPi85u88jZ/d5hJua6Fqzlqxddxn121oMBrUsutaspf3jpXQs/QR3Xi7Z8+Yw7kdXjNhCznhiFoka0gJPQUHS4jPtS2xJnryDD8S/w7Thb1hvWViBAA1PPEX70mW4c3PJnjuHyit/EJPqRlJQ0MRMIBwVGMdjSDiS4be3maipRfx+8g89iOzZuw2qq6eWRWDtetqXLKV9yTLC//3fuPNyqfxhCjqbnoyO4HJCMI7HkHD8O0xDLYvKH//Q1iB79XU2P7eQvIMPJHefvXod5g+sW791I/ysbLLnzaHiisvwVZRTcORY2Kx+RCYQpgzG8RgSjsvnxT91Cl0rV5M9Zzeyd5tF17r1NL+6iKZ/vUTufvuSd8D+hBsbbWezZBmurEzb2Vx+2aC0yFMS0+IxGOJL5s470fn5F2TP2Q0A/5TJ+L9zAaG6epoXvUHVr3+Lp6jIdjY/uLTPBZ9jhtEzgTAhGMdjSAqZM3eiedEb2y3q9JaWUHLaqRSfetKYGd4fNGnU4jHrpQxJwet0l0I1tb1eTzunA9jz72I5UgfjeAxJQUTI2nkmnZ99kWxTRg8a45FCmK6WYUhoOEyorp5gdTWhTTUEq6uJNDXhKSzEW1ZG5MyzCaz/Cm952YBLEzJn7kTr2++Sf0jfm6+nFSnmPGLBOB5Dr2gkQqiunlB1NcFNNfZrdQ3hhgbbyVSU46uoIHvObDyFBYQ3NxGqqbVXrT/2BKG6Olz+TLxlpfZRXua8L8NTVIi4XGTuOJ26hx7FCoZw+VJfsTQmzARCQzqjqmz6020Eq6pw5xfgqyzHW1FB1m6zyD/ycHtPnX726fEWFzH+up/Yks3NLYRqagnV1hKqraPj089t59Tehqe4GG95GeJy0bV6NVk7z+yzzLTBtHgM6YqIEG5oYPxPr4lpjZO4XHgKC/AUFpA5c8dtrlmBgN2aqq3DN34c3uL4ScakFGY43ZDOeIqLCTc3x21xpSsjg4wJ48mYMD4u5acqkkYtHjOqZdgOb0kx4br6gTMaRo5YR7RSzGkN2vGIiFtEPhaRfzjnh4nIRyKyRET+LSLTnfQcEXlWRF51ZDAQkQtFxBKR2VHlLR9gt3tDkvCUlhCqb0i2GWmG2F2tWI4UYigtnquAz6LObwfOUdW5wMPAz530c7H11K8CrozKXwX85/BNNSQKb3Ex4XrT4kk4psWzLSIyATgOuCcqWbFVBQHycTaEZqt2usW20yn/AewqIjvFYrAh/nhNiyc5pJHjGWxw+Y/AdUBuVNolwEIR6QRagH2d9IeARwA/cF5Ufgu4CfgP4ILBVGq005NTn5aPJzhxKoIgGRm4MjJwZfgR98C/U6P92UZLXb1rfI9I0SnBYJREu3WWPxSRg6Mu/Rg4VlXfE5FrgZuBSxxJ02P6KO5h4D9FZOpgjDPa6cmrT7EIVm2ia+Uqulatpmv1Wtz5efin70DmjB3wT98Bd07OiNQVC4msL651KUgCJhAOVzvduXYBW0Mqv1bVB5z0PbClkTOxN4u/ytm7uU8G0+LZHzhRRI7FbsXkicjzwExVfc/J8xjwr4EKUtWwiPwvtuCXYRQjLteWIe/8Qw5CLYvgho10rVxF63uLqX/0b7gLCsicvgP+GdPx7zANd052ss1ObRLT4pnPMLTTRaQIW5FiT2xLPxSRZ1V1s5PnUuBdbMdzNPDP/owYjJLo9dhKgTgtnmuwtXOqRWRHR1P5CLYNPPfHfLbvthlGOeJykTFxAhkTJ5B/6MFoJEKwagOdq1bT+s571D38GN6iQsK//BWR9nbc2cYJjUaGq52OLYnzkqo2AojIS8DRIrIIyFPVd5z0Bdj+ITbH04fxYRH5LvCkiFjAZuA7g7w3KCJ/xm7OGVIUcbvJmDyJjMmT4LBDtjgiVNnwP/9LyRmnkbXrzsk2M6UYJRMI+9JO7y+9qpf0fhmS41HVRcAi5/1QtNPnY7d0us//DPx5KHUbRjfdjshTWEDpeWdT//BjdHyyI0UnnzBqBfRGHbHPxSkRkcVR53ep6l1DLKMv7fShpveLmblsGHEyZ0xn/E+vRi1lw00307lqdbJNGv2MzMzlelXdM+oYqtOBvrXT+0uf0Et6vxjHY4gLLr+f0rNPp/jUk6h74CEan3kOKxRKtlmjm9Exj+dZ4Hyx2RdHOx14AThSRApFpBA4EnjBudYqIvs6I2LnA88MVIlZJGqIK1mzdmX8lMnUP/4km/54KxU/uNQEnvsgETGe4Wqnq2qjiPw38IFT1A3dgWbg+2wdTv8nAwSWwTgeQwJw5+RQdtH5bH5uIdW33UnF5Zfhzs5Ktlmjj1Gsne5cuw+4r5f0xcB2c4L6w3S1DAlBRCg84Vj8O86g+va7iHR0JtukUYUoiBXbkUoYx2NIGCJC0UnH4582lZo77sLqNM5nG8zqdIMhPogIRaeciG/SRKrvuAerqyvZJo0eRkdwOSEYx2NIOCJC8akn4xtXSfWd92IFAsk2aVQgGtuRShjHY0gK4nJRfNqpeMtKqbnLOB/AtHgMhkQgLhclZ3wbT1ERNXffjxUMJtuk5BFja8e0eAyGISAuFyVnnY47L4+ae+an9yRD0+IxGBKHuFyUnnMG7uwsau+dj4bDyTYpORjHYzAkFnG7KT33LCLNLXStWpNsc5JCOnW1zMxlw6hB3G4QwZWus5pTzHnEgnE8hlFFpLUVd24a7hGXgq2WWDCOxzBqUMsi0taOO3f7vZzTAuN4DIbEE2lrw5WVaXe50hHjeAyGxBNpacWTlzdwxjGIkF5drVgkjOeLyFpHwniJiMx10l0iskBE3haRXZ20g0VEReSEqPL+0UMux5DmRFpacOelYXynGzOc3is9JYwBrlXVuc6xxEk7EngPOAW4OiqvkTA29EukpRV3mrZ4zMzlXuhDwrgv+pIwXgo0i8gRQzXSkB5EWtJ0RKsb0+LZjm4J457bDd0oIstE5BYRyXDSXgAOwt679eYe+X/NViVCg2Eb7K5WmrZ4IK0cTywSxtcD1YAPuAtbHfQGVQ0DZ/ZWlqq+KSKIyAGDMS4VNbHTvb5Y6gqddTauzEzcWYOfQJgqz9aT3rTTU20XwVgYroTxg6p6rnM9ICL3YyuMDoYbsWM9Ay7IGROa2GlWXyx1bfzjrRQefwyZ03dISH1DJd7a6anWaomFAbtaqnq9qk5Q1SnYLZlXVfVcR9a0W+T9ZGD5YCpU1ReBQmDOsK02jEkirek7nA6JCS6LyNEi8oWIrBKRn/VyfbKIvOKEUBY58V1E5JCoEewlItIlIic713od4e6PWObxPCQipdgB5CXAZUO490YGob1jSB9U1Qynx7nFIyJu4DbgCOxR5g9E5FlV/TQq2x+wtdMfEJFDgd8C56nqa0D3lJkibPmbF6Puu1ZVnxisLbFIGB86nPuc82ehV+lTQ5qizg6EkpExQM6xSwKGxPcGVqnqGgAReRQ4CYh2PLsAP3bevwY83Us53wb+qaodwzXEbIthGBV0z+Gxe+5pSuyjWiUisjjquLRHDeOBr6POq5y0aJYC33LenwLkikhxjzxnAo/0SOtthLtPjOMxjArCLS1mDk/sjmcg7fTevHrPdtY1wEEi8jH2tJgNRA0EObHd3bCnzXRzPTAT2Asowh7h7hezVsswKoi0tOLOT+PAMgmJPVQBE6POJwAbozOo6kbgVAARyQG+parNUVlOB55S1VDUPZuct4Me4TYtHsOoILD+Kzz5+ck2I7nEfwLhB8AMEZkqIj7sLtOz0RlEpEREuv3C9WwvWXwWPbpZwxnhNo7HkHTaPvyY9qXLyD/s4GSbklTiLWHsTO69Arub9BnwuKquEJEbROREJ9vBwBci8iVQjj0CbdsnMgW7xfR6j6IfEpFPgE+AEuwVCv1iulqjEKuzk6516/EUFeEtLkI8Y/fP1LV6LQ1/f5rKyy8zLZ4ETCBU1YXAwh5pv4h6/wTQ67C4qq5j+2D0kEa4uxm7/9EpiNXVRcsbb9G86A28FWVEmluJNDXhLsjHW1qKt6wUb2kJ3tJSPGUleAoKEFfqNlpDtXXU3r+A0vPOxjeuMtnmJJcUXGEeC8bxjAKsYJCWN9+i+dXXydxpBuN+dAXeslIANBwm1NBIuK6OUF09wY2baF+6jFBtPVZHB56SYrzFxeBxg6rT33c6/QqghL7/A6rvvGfLOQqu7Gx8FWV4y8vxlpfhLS1J6M5/kfZ2qu+8l4JjjyJr5k4Jq3dUYxyPIRFYoRCtb71L8yuv4Z82hcorLsNXWbFNHvF48JWX4Ssv2/7+QIBQfQPh+gbUiiAIiHPAlvfu7Bzy9t/PGTYREHsUKVRTS+t77xOqqSXS1IynuMh2RBVl+ByHBGB1dGJ1Rh9d9msgAC4X4nbbTsvtJnzW2Wxe+C9w0rrTxeWyu4xuN+J20fLm22TPmUXefvvG+VNOPlZXF7jduLzefvOZFo8hrqgqLW++RdNLr5IxaQLll11CxvhxQy7HlZFBxvhxA97ryvSTNWuXfvNYoRDh2jqCNbWEamroWL6C0KuLbLmZTD+uzMytR1YmnrJSXD4falkQiaARC6wIiAvEhYbDaCCAhiOoFYGIhUYiaDgCVgT/tKkUHnf0kJ95tKOWRaimlsC69XStXU9g3TrCm5vAsnDlZOMtKcZTUtLHzYm1NZkYx5NANBKh9b0PCO67Px2ffkb5dy8iY2LitsfoD5fXi2/8OHzDcIDRePJyKTzmyBGyavRjdXbStf4rAmvXEVi3nsD6r3FlZ5ExZTL+KZPJO3B/uxUrQnhzE+H6BkINDbS9+/52ZZkWj2FE0UiEtsUf0fTCS3hKSvAefyIV37sk2WYZhkGoto6uNWvpWreewFq7NZMxcTwZU6aQ+839KT13Up8zsL3F9ihlJjNoeKzHwFGabYthHE+cCbe0UP2X23Hn5VJ6zpn4d5iGy+dLtlmGYaDhMBv+909k7bqz3ZrZ/xv4xlWOXFDeOB7DSOHKzMSVk4Nv3Dgypk1NtjmGGLCCQTJnTKfs/HNGvGwhvXYgTN1JICmCy+ul4tLv0LVmLU0LXxj4BsOoRUMhfBNii4H1h6jGdKQSxvEkAFdmJhXf/y7tS5fR9MpryTbHMEw0GMI3fruJuyNU+AgcKYRxPAnCnZNDxQ8upfXfbxPp7Ey2OYZhYAWDZEyIk+PB6GoZ4oSnoIDseXPQ0ID73BtGGeGWFgDcBXFcT5ZGLR4TXE4wkfaOlF5flbaEI2BLM8WtilRrtcRCLNrpDzm71S8XkftExOukG+30frA6OsBtHE+q4crNsWdoW3EcekqjFk8s2ukPYW93uBuQCXTPiDPa6f1gdZgWTyri8nrB5cKKV3wuxvhOqrWWhq2drqoL1QF4H3sbRTDa6f0Sae8E43hSEnG5iLS2xq+CNGrxDDbG062dvt1ccKeLdR52iwjs3c0eBM4Heu5y/2vneGkwlaaiNO1ABK67jorMLCNhnIL1lYqL64oKcY2ABE9PCWN7AmGKeY8YiEU7vZv/A95Q1Tdhy/aKI6KdPiakaXuw/j9/wR/vvYdbajcOnHmEGDMyv0mqT1Vpfec9rt9nP26prop5IW1fpFp3KRZi0k4XkV8CpcD3hlDnoLXTxyKewgJ7awhDShBpbaXukceJNLfgPeFkfO44BZdTsLsUC7Fop18CHAWcpaqD/muku3a6p7gIjaSlz005OpZ/yoabbsZXWcm4H/8Qlze+s0/ivdn7aCKWKOcd2LvQv+MItf9ioBuiuJGtwei0wltWtkWu1zA6sQIB6h9/koYnn6L0gnMpOuHYxGy4n4Dgsogc7UyDWSUiP+vl+mQRecVRBV3kDCx1X4s43/UlIvJsVPpUEXlPRFaKyGOOdE6/xKKdPuh7jXb6VvIO2J9Ieweh5ka8RUXJNsfQg8BXX1P314fJmDyJ8df9BFdmZsLqjneMR0TcwG3AEdjTWz4QkWdVNVo7/Q/AAlV9QEQOBX6LPXgE0Kmqc3sp+n+AW1T1URG5A7gYuL0/W8y4boLxFOTjzslh8z/+lWxTDFFoJELTiy9Tfec9FBxzFKXnnpVQp7Nlk/5YjoHZG1ilqmtUNQg8CpzUI88uwCvO+9d6ub4NjojfoWyVxHkAW9SvX8ySiSTgzsula9Uqghs3GVmXBKPhMKG6ekLVNQSrawjV1BDcVE24oQH/DtMYf+2P8RQUJMW2EWjxlIjI4qjzu3rop48Hvo46rwL26VHGUuBbwJ+wJwHnikixqjYAfqf8MPA7VX0aKAaanNHs7jIHXElrHE8SEBEypk4hWFNjHE8cCVZXE9xYTai62nYy1TWEGxvxFBY6ahrlZM3alfzDD8VbVobL178KRNyJ3fHUq+qe/VzvLbzRs9ZrgFtF5ELgDWADW0egJ6nqRhGZBrzqqIe2DKLM7TCOJ0m4MjLQLhNkjhedK1dRe98C/DN2wFdRTvbc2fgqKvCUlgwoM5MMRDUREwirsCWIu5kAbDOhTFU3AqcCiEgO8C1VbY66hqquEZFFwDzgSaBARDxOq2e7MnvDOJ4k4fL7sYzjiRudn31O3oH7U3jMUck2ZdAkYALhB8AMEZmK3ZI5Ezh7GxtESoBGZ4rM9cB9Tnoh0KGqASfP/sBNqqoi8hrwbeyY0QXAMwMZYoLLSUI8bqxgMNlmjFk6v1hJ5o47JtuMoRHn4XSnRXIF9rKmz4DHVXWFiNwgIic62Q4GvhCRL7Gny9zopO8MLBaRpdhB599FjYb9FPiJiKzCjvncO5AtpsWTJDq/WEnRSccn24wxSdvijwg3NZMxeeLAmUcRiVgyoaoLgYU90n4R9f4Jto5QRed5G3snit7KXIM9YjZojONJAlY4TKSlBf/0HZJtypij+bU3aF70OpWXfy8xk/5GCgXMIlFDPLHaO8ieO8fsyzOCqGWx+bmFdKz4lHFXXYGnqDDZJg2d9PE7xvEkGlXF6ugge495yTZlzKCRCPWPPE6orp7Kqy7HnZ2dbJOGRTqtTjc/uQkmWLUBgIzJk5JsydhAVam5+z4iHR1UXP69lHU6QCJmLo8ajONJMJ2ff4Er0x/XTcPTiXBDA+7cXMovvjDlpaHTaetT09VKMMGNmxCvD7qSbUnqE/i6Cmv8ZErO+PbI6ZcnCdH02oHQtHgSTHDjJiTZU/PHCE0vv4onNze1Rq/6w4rxSCHGyF8sNbA6Owk3NIydL0oSCdbU0rVqNa6cbGhvSrY5I0Kq6Z/HgvkGJIhwSws1d9xD7n7fMPGdEaD5ldfIO2D/sfNZmq1PDSNNqL6eTX+6jazdZlF0yokD32Dol/DmJjo+WU7eAd9MtikjSIwjWinWWjItnjgTqNpAzZ33UnDU4eR9c79kmzMmaH7tdXL22Rt3dlayTRlRUm1kKhZMiyfONDzxFIXHHmWczggRaW+n7YPF5B9yYLJNGXnSqMUzoOMREb+IvC8iS0VkhYj8Pyf9UBH5yNFOf0BEPE660U6PItzQSObMnZJtxtghYoEIVkecpISThRqViZ4EgENVdQ4wFzhaRPbD3lv1TFWdBazH3ocDjHb6FjQcJtLejjs/L9mmjBncebkUHn8sdQ8+gkbGmD6ZafFsxZFHb3NOvc4RAQKq+qWT/hL2Pq1gtNO3EN7chKcg3ywGHWFyv7EP7txcml54OdmmjChiaUxHKjGo4LIji/EhMB1bHuN9wCsie6rqYuzdx7o3PzHa6Q5WQSnhX/0/fGWlCamvL8aidrr+138RrK6hzONNyWfrqZ0OpFyrJRYG5XhUNQLMFZEC4ClgV+xtE28RkQzgRZwNoY12+lZa332frtVrKD1n249jLOqLJ6OutiUf8zO3iz8HOxOyUXtcn01JudnHsTCkPoCqNmEL8x2tqu+o6gGqujf2bvQrB1lMt3b62EeVSEsLmka/ZInEnZONhiOENm1KtikxI6i94XsMRyoxmFGtUqelg4hkAocDn4tImZOWgb3n6h2DqTCdtNNz9tqDcFMzHUuWJduUMYVaFptfeJm6Bx/FW1oydrYYMcHlbagEXhORZdi71L+kqv8ArhWRz4BlwHOq+uoQ6k0L7XTxeCg549s0PPUMkbE2/JskIu0d1Nx9H52ffc64a67ClZGRbJNGjjRyPAPGeFR1GbZ+Ts/0a4FrB1NJOmun+6dNJWvWrtTeN5/seXPJmDQRLR9QaNHQC4Gvvqb2/gVkzd6NohOPS/mtMLYhQTEeETkaWyXUDdyjqr/rcX0ytqRNKdAInKuqVSIyF1sPPQ97VPtGVX3MuWc+cBDQ7BRzoaou6c8OM86bAIpOOp7sObMJrFtP3YOPENywkbqHH0u2WSlFpKOT6tvvouik4yk+5cSx5XQc4h3jcUanbwOOwdZIP0tEdumR7Q/AAlWdDdwA/NZJ7wDOV9VdgaOBP3aHYByuVdW5ztGv0wGzVishuDIyyDtg/y3nvrJK2j9aQsmZp5k5PoOkY/kK/NOmkT13DIcG499d2htY5cjRICKPAicBn0bl2QX4sfP+NeBp27Qtc/ZwZIxrsVtFw9qTxPzXJwFxuXBlZRHevDnZpqQM7UuWkj13drLNiB+qYFmxHVAiIoujjp7z6MYDX0edVzlp0Sxl62TgU4BcESmOziAiewM+YHVU8o0iskxEuqfY9ItxPEnCW1ZKqLYu2WakBFZnJ12r1pA1q2evYIwR+w6E9aq6Z9RxV48aeour9mxmXQMcJCIfY8dtNuDM0QMQkUrgr8BFjswx2FLHM4G9gCLsUe5+MY4nSXjLSgnVGMczGNqXf0rmjOm4MjOTbUpcScA8niq2rjAAe2R5Y3QGVd2oqqeq6jyc+Xaq2gwgInnA88DPVfXdqHs2OUurAsD9DEJV1DieJOEtLyNUW5tsM1KCjrHezeom/sPpHwAzRGSqiPiwVxg8G51BREpEpNsvXI89woWT/ynswPPfetxT6bwKcDKwfCBDTHA5Sbizs+lqaxs4YxpiBYN0rV5L15cr6fxyJZGWFkrPPSvZZsUXJe4SxqoaFpErsNdTuoH7VHWFiNwALHamuRwM/FZEFHtFwuXO7acDBwLFInKhk9Y9bP6QiJRid+WWAJcNZItxPEkiuKkaX2VFss0YNQSqNtCx4jO6vvySwFdV+CaOJ3PHGRR/62QyJk8ak8Pn25KYSYCquhBY2CPtF1HvnwCe6OW+B7EXf/dW5gr95JgAAA9cSURBVKFDtcM4niQRrNpA3gFmV0KA9o+X0vDk02TvMY/8ww7Bv8O0sTUjebCk2OzjWDCOJwkoEKiqwjdhzK8aGZDOz7+k/omnqPjBpWSMH5dsc5KLcTyGuBKJADLknQnVssbUhMOudeupXfAQ5RdfYJxOAmI8ownjeJKAFQyRMXH8oDWhQg0N1D/2JIHVa/CWleKtqMBbUY6vohxfZQWe4qKUi4EEq2uovWc+pWefgX+Hack2ZxSgYI2xrVz7wTieJKDBIJ7CQjqWryCw/mu61n9FqLoG3/hKMnfckcydZuCtrABVWl7/N00vvUz+YYdQftF5hOrqCVbXEKquofW99wlV1xJpacFbWoK3spKik4/Hkze693gONTZSffvdFJ50/NifFDhYTIvHEG80EqH94yWEGhrxT55E/oH7460oJ1i1gc4vVtLy5ltoMIgrKxN3bg7jfvRDvM72qRmTJpIxaeI25VnBIKGaWlpef5OmF1+m5NunJuOxBkWkrY2a2+8m/5ADyd1rj2SbMyQUCDc3Y7V3YAUCaCCIFQg47wNR74NYXQE0GMDqCmAFA2hXACsQxJXhI3PnmX1UYByPIY54igqZ9JsbtovXeEtKtiyCDNU3EG5owD9j+oBxHZfPR8bECRSddDxVv7mJgsMPw1OQH7Od4c1NaCgEbhfidiMuN7jdiMdtd+1cLsTlQi0LVFHLQlWxurrsNMtef6Rqv9dQiNoFD5E1dzb5B48uXSyrq4twcwuR5mb7tanZee+8NrUQ/K//YsNNt+DOycaV4UMy/M5rBi7nkAwf7sJCXH77vZ229XqkvZ2OTz/r3QjjeAzxRGBAZ+ItKcZbUtxvnp64c3PJ2Xsvml95jeL/3965R1ldXXf8871z5/24AzMwgyAIQgQh4EqMITQS6ysYH6hRY/OoCdEVW7UkMVZdzaP1EQXSZtEmWq200rXyaCpVgwSRGK0VBwMmSkAE8QEMFshyBspw53nv7h/nXLiD87hz5wEzcz5r3TW/3++e39nnwP3t3z7n7LP3Zy/Pun2tdXXUr1xN4xvbiBQVQTKBJdyHRBJLtGGJpJ8k90QiINGyeDG7vns3iggUgYhcX/1x8ayZjPjMvKzb1hsskaBpx9s07thBov7AUUVz4CAkE+TEYkRjMfe3vIxoZQUFp04ipzxGTqyMvHFjmXDv3/aqDbmMomDiKRx85tljWxcUT2DwEjvvHPbct4TY+ecS7eGqWbKpiQO/fo5D62oom/snVF57VZf+NOZd9dOVaH71WE5ZfG/W7e9rki0tNL6xnfimPxDfspVoxUiKpk2lYPKpXtGUkVMeI1JY2O1kf6aLAVlhpHaYDwuC4hliRMvKKPnYmc7quXJ+RvdYMknDyxuo/9UaCk+bwtjbv0m0vLzb+yRBfz6MWZI4HCf++lbim/5A4/Yd5J88juKZMxhx8UVER3Tfr+NGsHiOIqkAt2cj35d/zMy+5zeE3QNcjQuF+KCZ/aPfYPYoLgfXDX4vyDm4oEKXmdlKX+9TwA98WNRAHxI770957++XEq2sIDb3k12WTcQb2ffwMoiIqhu+8oGJ68FEc+0e6leuoundXRROmUzRzOlUXns1OcXFx7tpmREUTztSKYwbJOUCL0paDUzDbbGfambJVNYJjqYwvg0XNnGBv55KYbyyLzsQ+CDRWBljFv4lex/8F5INDdiC6zssl4jH2fvAwxRMmsjIKy7r36FEP2JtbdSvWcuhl9Yz4uKLGL3gukG45cLCcno65pJCHZvC2IC/AD6fCgZkZqkYD12lMM6VdIGZZZRJNJA9uRUVnLTwJvY+tIy2+nqXwz3tzZ84HGfvAw9RMPlURl5+6aBVOsmWFvYs+SG5o0Yx9q9v7fG81gmDMfRywXdBVimMzexlSacCn5N0BfBH4K/M7E36MIVxoHfklJYy5uYbIZlk9133EcnLJbe6mrwx1TTteIuCD01h5PxLBqXSSba0UL/qaVrnX075vAsoPmPWoOxHO8JQqz3HpjCWNAM359NkZmdKuhIXMOjsvkxhPBhzYp+o8pYufxQSCay11fnmSERKSvo8x1B/9M2SSay1zbW9rRVrbSOZEyFyzbWMra7iuycNzD6vfs2dnoq5PEzo0aqWmR2Q9DwuvUUtsMJ/9Tgu5GEmpFIYt3VXcCjm+x5oec27a/nWxMn8zQM/RnK+Ns6RSEii8PRplJx1Zp9ZC73tW9O7O2neuYvWfftp3buPln37sLYEeVWjya0a7faoVY12+9UqK/hma+vQ+Z0Ei+coPrJYq1c6qRTGi3BpL87FWTqfArZ3XstRzOwZSXcDw3w7cv9jZrz/X0+iO+6kbM5s53eTTDmqGdaW4OCzzxHfvIXKz11NTsnxW/1p2befuidX0rp3H4VTTyO3uoriWR8mt6qKnFjZ4B9GZYAFi6cdY4Dlfp4nAvzCzJ6S9CIu5OE3cJPPHS+ddMy9wJM9bm2gRzRu207y8GFyYmUUzZjeYZniM2ZSv2o1exb/A5Wfv4aiqacNaBsT8TgHnl5Lw8ZXKD//XKoWXIeiw9G9LHgut6OLFMYHgIszETKcUxgfL8yMA796mvJ5F3b5D61olJHzL6Vw2lT++JP/cI52l15MJC/3SD3W0kIy3kji8GG3QTIexzDyTx5HtKIiK2vEEgkOrauhfs1aimfNZNydt5FTWpplb4cAYXd6YCjQuGUryZa2jLMzFH7IeSy//4sV1H5/EZGCApLxOMnDcTcRXVxETlERkeIiIsXFkEhS98RKrLWV/PHjyZ8wnvxTxmOjux9Bx1/fSt0TK8mJxRhz043knTSmt90dGlgYagUGMWZG/eo1jLjowh5FLMwpKmLUdV+kpXaPy3ZaXEykqOiI9dMRbQcP0rxzN807d3Hw2edpmTKN3ffdT/6Ek8mfMJ6cslIUjbrd7Ab/98KLtNXVOStr+rRhMXeTCQbYAFg8kuYBS3H+do+Y2f3HfD8BN287CqgDvmhmtf6764Bv+6L3mNlyf/2juN0KhbhA8gu9/1+nBMUzBGnY+DuQKJo5o8f3SiL/5MyXjKOxGNGZMYq9rLzqsVTd8GWa391F867dNL/9jtvZ3uZ2txdNP52ys+cMuoiJ/Y5ZvzsQ+nnaHwMX4FalN0j6pZml507/AS531nJJ5+J2H3xJ0kjge8CZOD35ir+3HngQ57O3Hqd45gGru2pLUDxDjERDA3VPrKT6a189LtaEgLzqavKqqymd3W1CyUA6/T/UOgvYYWZvA0j6OTAfSFc8pwPf8MfP4VavAT4NrDWzOn/vWmCed68pM7Maf/3fcUn9ulQ86sYiOm74hGKBwFBlp5mdkjqR9DRQ2cs6C4CmtPOH0/OnS7oKmGdm1/vzLwEfN7Ob08r8FHjZzJZ6x+AVvl1fAQrM7B5f7jtAI27R6H4zO99fPxu43cwu6aqhJ6zFY2Zh8B8YNpjZQERH6+iZOvYF/y3gRz5b6AvAHpyzb2f3ZlLnBzhhFU8gEOhzanERJVKMA95LL2Bm7wFXAkgqAT5rZgcl1eLSG6ff+7yvc9wx19vV2RFDJ0lTIBDojg3AFEkTJeXh9lT+Mr2ApEofUwvgTtwKF7jN3xdKGiFpBC78zRoz+1/gkKTZPkbXn5OBc3BQPIHAMMFv4L4Zp0S24nYhbJF0l6TLfLFzgG2StgNVuF0G+Enlu3HKawNwV2qiGRci5xFgB/AW3Uwswwk8uRwIBIYuQ87ikfSvkvZL2px2bYmkNyRtkvS4D++R/t1GSZ/y549Lujzt+22Svp12vsLP9ncmf6GkzZK2SPq6v3aSpN9IelJSiaRySe970xRJn5Bkksb585ikujSTt6v+lkt6zPdvq6+rz+VJKpD0W0mv+b79nb8+XVKNpOWSIpJmSXo17b4/kxSXi16JpA9L2tRdv7poxzz/f7JD0h0dtaEXdWfUR3/tUUnvSHrVf17KVu6wxMyG1AeYC3wE2Jx27UIg6o8XAYv88VRgCVCEMzvBhWxd7I8rcAHQVqXV9R5Q3YnsGcBmX18U+DUwBbgfmA5cCtzoy24BTvfHtwK/A67x558GVmfY3+XA9f44DyjvD3m41YsSf5yLC287G1iG83K9Bec4FgHqgVJf9p+8rLP8+ddw8bmz+b/NwZnyk3xfX8P5nbRrQy9+Oxn10X//KHDV8f69D9bPkLN4zOwFnKt3+rVnzI1vwXlXpmbhU2Fa05cF1wFz/PEc4ClglBwTgUYz29uJ+GnAejOLe3n/DVxBx+Fgj5Xzw2POu32DSirDKdplvp8t5jbv9rk8c3QUAjfH/03ihu5J3BzAx33Zj+K8ZXvUt0444gBnZi1AygGuXRuyrDvjPmZbf+AoQ07xZMAC/OSXmW3BWScv4ty+wVk4M/ys/xygBtiGUypzcA9wZ2wG5kqqkFQEfAa3fPkj4CHgRlxYWHAPX+phnAT8J84dnQzkpJiECzv7b5J+L+kRScX9JU9Sjh9G7cd5sb6M2/ezCvgE8Ey6LN+WJG7ZNV3xZNK3jhgL7E47r/XXOmpDVvSgjwBL0oZaP+mN3OHGsPLjkZSKfHjkR2Jm7aJQmlmzpC244dpsYDHuQZ2DCw/S6dvazLZKWoSLKd2AGwq0mdlOnGWSzjrgDm9FvWtmTd6qKsFZCb/NoEtR385bzMXBXgrcYWbf6Q951kEIXDP7PUetm3RZtwL/A2wws7ckTZYLKldi3mU/Czp0VuukDVnRgz4C3GZmj/WF3OHGsLF45HbWXgJ8wfwgvQtewj24peY2wa3HKZ5u39ZmtszMPmJmc3FDvjc7KfcmMAI3D1PjL7+Cc01/J83k74paoNa/lQEewymi/pKXqusAzorpzNt2PfAx4JNpsmpxfiO9mYTt1gGur8igj4FeMCwUj1wogNtxCQXjGdyyDjcJ+po/34SzfsbjJmm7kjXa/x2P8wD9WRfFa4CFHH04a4Cvk+HD6eeadktKhQ08j/Yb/vpMnqRR3gpAR0PgvtFJuw7hhkRfzkZWF3TrANcbetLHQO8YcopH0s9wP/LTJNVK+ipuzqMUWOvH4//cTTUv4YZXNXDE8Wo/sNFPnnbFCkmv4xIX3uQtps5Yh3uDb/TnNV5uTx7OW3AhaDcBZwDf7yd5Y4DnvJwNuPmPp7qRlW9mqTmZbPrWDuvEAS7b+jqgp31Mn+N51SvDQAYEB8JAIDDgDDmLJxAInPgExRMIBAacoHgCgcCAExRPIBAYcILiCQQCA05QPIFAYMAJiicQCAw4/w/WqZlLouYoGwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fieldset = FieldSet.from_data(data, dims, mesh='spherical')\n", "fieldset.U.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, printing the velocites directly shows something perhaps surprising. Here, we use the square-bracket field-interpolation notation to print the field value at (5W, 40N, 0m depth) at time 0." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1.1747725785927632e-05, 8.999280057595393e-06, 20.0)\n" ] } ], "source": [ "print((fieldset.U[0, 0, 40, -5], fieldset.V[0, 0, 40, -5], fieldset.temp[0, 0, 40, -5]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While the temperature field indeed is 20C, as we defined, these printed velocities are much smaller. \n", "\n", "This is because Parcels converts under the hood from m/s to degrees/s. This conversion is done with a UnitConverter object, which is stored in the .units attribute of each Field. Below, we print these" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U: \n", "V: \n", "temp: \n" ] } ], "source": [ "for fld in [fieldset.U, fieldset.V, fieldset.temp]:\n", " print('%s: %s' % (fld.name, fld.units))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the U field has a GeographicPolar UnitConverter object, the V field has a Geographic Unitconverter and the temp field has a UnitConverter object. \n", "\n", "Indeed, if we multiply the value of the V field with 1852 * 60 (the number of meters in 1 degree of latitude), we get the expected 1 m/s." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0000000000000002\n" ] } ], "source": [ "print(fieldset.V[0, 0, 40, -5]* 1852*60)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### UnitConverters for mesh='flat'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If longitudes and latitudes are given in meters, rather than degrees, simply add mesh='flat' when creating the FieldSet object." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEjCAYAAAChAAFiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXGWd7/HPl7B5JcgSxLBIwsjo4BYxIg6OF1ARuAgojoILiN5hdEDQubgwepHB4YU6rowOGJVNEGRwkFyNIoMEHBUEZAmLSFiUhiAigigaTfK9f5yn4aRTXXU6larqTn3fvM6r6+y/qi76l+c5zyLbREREDKt1Bh1ARETEICURRkTEUEsijIiIoZZEGBERQy2JMCIihloSYUREDLUkwpgyJD1T0nWSHpV01KDjmawk7SZpZNBxREwVSYTRM5Is6Rljth0v6ezVvOT7gIW2p9s+WdIZkv6l+0hXik+SPibp12X5uCS1OX6OpGslPVZ+zpkM14qI5pIIYyrZDrh5TV1M0rotNh8OHAA8H3gesC/w9+Ocvz5wEXA2sClwJnBR2T6wa/VLSc75GxJTXr7EMWlI+gtJ3yslngclnSNpk7Lve8DuwOck/U7S4cCbgPeV9f9XjttK0tcl/UrSXfUq1FIavUDS2ZJ+C7y1RRiHAp+0PWL7XuCT4xwHsBuwLvAZ20ttnwwI2GPA1xp9v/9H0gOSlkg6rLb9KZLOKp/RzyV9aDShjS2xS5pVSvbrlvWFkk6U9APgMWB7SW+VdGepsr5L0pvaxRUx2SQRxmQi4CRgK+CvgG2B4wFs7wF8HzjS9ka25wHnAB8v668uf8z/H3ADsDXwcuDdkl5Vu8f+wAXAJuX8sZ5dzh91Q9lWBSjdKOmNtWNv9MrjFN5YO75v12rhacBTqD6HtwOfl7Rp2fdvZd/2wP8EDgEOa3WRcbyFqoQ6HfgVcDKwt+3pwF8D10/gWhED16pqKGIgbC8GFpfVX0n6FPDhCVziRcAWtk8o63dK+iJwEHBx2fYj298or//Q4hobAY/U1h8BNpIkV57X5tjR46f3+1ot3sefgRNsLwMWSPod8ExJVwNvAF5g+1HgUUmfpEpuX25xnVbOsH0zgKRlwArgOZJ+YXsJsKThdSImhZQIo5eWA+uN2bYe1R/pVUh6qqTzJN1bqi7PBmZM4H7bAVtJenh0Af4J2LJ2zD0drvE7YOPa+sbA78ZJNmOPHT3+0QFfC+DXJQmOeowqmc4A1gd+Xtv3c6qSY1OPf4a2f0+VWN8BLJH0LUnPmsC1IgYuiTB66RfArDHbZrPyH+G6kwADz7O9MfBmqurS8YxNAvcAd9nepLZMt71Pm3PGupmqQcqo5zN+A52bgeeNab35vNrxg7pWOw9S/UNku9q2pwP3lte/B/5Hbd/TWlxjpc/Q9sW2XwnMBH4KfHE14ooYmCTC6KWvAR+StI2kdSS9Ang11TO6VqZTlXwelrQ18N4O1/8l1XOuUT8Gfivp/ZKeJGmapOdIetEEYj4L+EdJW0vaCvg/wBnjHLuQqtR7lKQNJB1Ztn9vwNcal+3lwPnAiZKmS9oO+Eeq0jdUz/deJunpkp4CHNvuepK2lLSfpCcDS6l+f8snGlfEICURRi+dAPwQ+G/gN8DHgTfZvmmc4/8Z2Inq+de3gP/scP0vAzuWatBvlD/yrwbmAHdRlX6+RNUwpKkvUDW4WQTcVOL4wuhOSTePtoq0/SeqLg2HAA8DbwMOKNv7eq0JehdVye9Oqt/NV4HTShyXUP0D5kbgWuCbHa61DlVSvg94iKrxzT+sZlwRA6FMzBsREcMsJcKIiBhqPU2EkjYpHZh/KulWSS+RtJmkSyTdXn5u2vlKERGxtpF0Whn0oeXjEknPkvQjSUslHTNm316SbpO0WNIHattnS7qq5Jiv6YnRmcbV6xLhZ4Hv2H4WVSu3W4EPAJfa3gG4tKxHRMTwOQPYq83+h4CjgE/UN0qaBnwe2BvYEThY0o5l98eAT5cc8xuqASXa6lkilLQx8DJKJ13bf7L9MNXIHmeWw86kaiAQERFDxvYVVMluvP0P2L6aVfse7wwstn1naVB2HrB/6X60B0+0TG+UY3o5ssz2VMMvnS7p+VQt0I4GtiyjT2B7iaSntjpZ1ViShwNo/fVfuN6WLQ+LiIg14E/3jDxoe4vR9Vft/mT/+qHuesJce+PSm4E/1jbNK8MjdmtrVh4cYwR4MbA58HBtMIkRGgwW0ctEuC5VU/h32b5K0meZQDVo+bDmAWzw9G299THv7k2UERHBXUcfs9JAF79+aDk/vvjpXV1z2szb/2h7blcXaa3VQBtus72tXj4jHAFGbF9V1i+gSoy/lDQToPx8oIcxRETEajCwosv/emiEalD+UdtQ9WV9ENhET0yxNrq9rZ4lQtv3A/dIembZ9HLgFmA+1ZQylJ8X9SqGiIhYXWa5V3S19NDVwA6lhej6VAPrzy9j714GvK4c1yjH9Hr2iXcB55RA76Sa6mUd4HxJb6cai/JvexxDRERMkIFlPR4tT9K5VHNxzpA0QjXbzHoAtk+V9DTgGqpB5ldIejewo+3flmEILwamAaeNzogCvB84T9K/ANfRYFaVniZC29cDreqHX97L+0ZERHeMWd7jkcdsH9xh//1U1Zut9i0AFrTYfidVq9LGMh9hRES0tKJzO5O1QhJhRESswsDyJMKIiBhmKRFGRMTQMvT8GeFkkUQYEREt9bQDxCSSRBgREaswzjPCiIgYYoblw5EHkwgjImJVRvy55dCda58kwoiIWIWBFSkRRkTEMFueEmFERAyrqkN9EmFERAyxFU4ijIiIIZUSYUREDDUjlvd07vbJI4kwIiJaStVoREQMrVSNRkTEUDPizx6OFDEc7zIiIiYsJcKIiBhatljuNJaJiIghtmJISoTDke4jImJCqsYy63S1dCLpNEkPSLppnP2SdLKkxZJulLRT2b67pOtryx8lHVD2nSHprtq+OZ3iSIkwIiJa6EvV6BnA54Czxtm/N7BDWV4MnAK82PZlwBwASZsBi4Hv1s57r+0LmgaRRBgREaswsKLHlYa2r5A0q80h+wNn2TZwpaRNJM20vaR2zOuAb9t+bHXjSNVoRES0tNzqagFmSLqmthw+wRC2Bu6prY+UbXUHAeeO2XZiqUr9tKQNOt0kJcKIiFjFGhpi7UHbc7s4v1VrncdnSZQ0E3gucHFt/7HA/cD6wDzg/cAJ7W6SRBgREaswTIYO9SPAtrX1bYD7auuvBy60/efRDbVq06WSTgeO6XSTVI1GRMQqTHfVosvXzDil84FDSuvRXYBHxjwfPJgx1aKllIgkAQcALVuk1g083UdExOTU68Yyks4FdqN6ljgCfBhYD8D2qcACYB+qVqGPAYfVzp1FVVq8fMxlz5G0BVW16vXAOzrFkUQYERGrsOl59wnbB3fYb+CIcfbdzaoNZ7C9x0TjSCKMiIgWNDQjyyQRRkTEKkzvS4STRRJhRES0lBnqIyJiaBllhvqIiBhuKRFGRMTQqmaonzboMPoiiTAiIlZhYEUay0RExDBbnu4T3ZN0N/AosBxYZnuupOOBvwN+VQ77J9sLehlHRERMjK2UCNeg3W0/OGbbp21/og/3joiI1ZR+hBERMbSqiXmHo2q01+newHclXTtmQsYjy6SJp0natNWJkg4fncxx+e9+3+MwIyJiZWK51+lqmSp6HemutncC9gaOkPQy4BTgL4A5wBLgk61OtD3P9lzbc6dt9OQehxkREXVVq1F1tUwVPa0atX1f+fmApAuBnW1fMbpf0heBb/YyhoiIWD3pUN8lSU8G1rH9aHm9J3CCpJm1iRVfQ4NJEyMior+MWJYO9V3bEriwmiSYdYGv2v6OpK9ImkNV8r4b+PsexhAREauhmo9w6lRvdqNnidD2ncDzW2x/S6/uGRERa85Ues7XjXSfiIiIVVSzT+QZYUREDLEMsRYREUNrtPvEMBiOcm9ERExQVTXazdLxDtWgKg9Iatl7QJWTJS0ug7DsVNu3XNL1ZZlf2z5b0lWSbpf0NUnrd4ojiTAiIlpagbpaGjgD2KvN/r2BHcpyONWALKP+YHtOWfarbf8Y1XjWOwC/Ad7eKYgkwoiIWMVo94luls738BXAQ20O2R84y5UrgU0kzRzvYFX99fYALiibzgQO6BRHnhFGRMQqjFi2ousO9TMkXVNbn2d73gTO3xq4p7Y+UrYtATYs114GfNT2N4DNgYdtLxtzfFttE6GkGxsE+ivbL29wXERETCFrYPaJB23P7eL8VgG4/Hy67fskbQ98T9Ii4Ldtjh9XpxLhNGCfDkHOb7M/IiKmoEnSanQE2La2vg0wOob16M87JS0EXgB8nar6dN1SKnz8+HY6JcK/t/3zdgdI+odON4mIiKlnEnSon081bd95wIuBR2wvKdP3PWZ7qaQZwK7Ax21b0mXA64DzgEOBizrdpG0itP3fnS7Q5JiIiJhi+jCVkqRzgd2oniWOAB8G1gOwfSqwgKpWcjHwGHBYOfWvgC9IWkHV6POjtm8p+94PnCfpX4DrgC93iqNRYxlJ+wIfAbYr56iK0xs3OT8iIqaWfsxQb/vgDvsNHNFi+w+B545zzp3AzhOJo2mr0c8ArwUWlcAiImItNwmeEfZF00R4D3BTkmBExHCYJI1l+qJpInwfsEDS5cDS0Y22P9WTqCIiYuCSCFd2IvA7YEOg47htERExtVUz1A+81WhfNE2Em9nes6eRRETE5OHhKRE2Tff/JSmJMCJiSIw+I+xmmSqalgiPAN4naSnwZ9J9IiJirTeVklk3GiVC29N7HUhEREweZmqV6rrRtmpU0tM6XaDJMRERMfXY6mqZKjo9I1zQ4BpNjomIiCmmDxPzTgqdqkafL6nVtBajROtpLyIiYgrzELUa7TTodtezMkZExNQ0lao3u5EZ6iMiogWxfEU61EdExJDKWKMRETHcXD0nHAaNy72SXirpsPJ6C0mzexdWREQMWlqN1kj6MDAXeCZwOtUMwmcDu/YutIiIGBSTxjJjvQZ4AfATANv3ScpoMxERa63hGVmmaSL8k21LMoCkJ/cwpoiImASG5Rlh00R4vqQvAJtI+jvgbcAXexdWREQM2rBUjTZqLGP7E8AFwNepnhMeZ/vfehlYREQMjt37sUYlnSbpAUk3jbNfkk6WtFjSjZJ2KtvnSPqRpJvL9jfUzjlD0l2Sri/LnE5xNG0sMxv4vu1LyvqTJM2yfXeT8yMiYupZvqLnJcIzgM8BZ42zf29gh7K8GDil/HwMOMT27ZK2Aq6VdLHth8t577V9QdMgmnaf+A9gRW19edkWERFrqV6XCG1fATzU5pD9gbNcuZLq8dxM2z+zfXu5xn3AA8AWq/s+mybCdW3/aXSlvF5/dW8aERGTm+kuCZZEOEPSNbXl8AmGsTVwT219pGx7nKSdqfLRHbXNJ5Yq009L2qDTTZo2lvmVpP1szy833h94sOG5ERExBa2BRqMP2p7bxfmtipWPhyVpJvAV4FDbo7WWxwL3UyXHecD7gRPa3aRpInwHcI6kz5XA7gEO6XSSpLuBR6mqUpfZnitpM+BrwCzgbuD1tn/TMI6IiOgHT4pWoyPAtrX1bYD7ACRtDHwL+FCpNgXA9pLycqmk04FjOt2kaavRO2zvAuwI7Gj7r20vbvQ2YHfbc2r/KvgAcKntHYBLy3pEREw27nLp3nzgkNJ6dBfgEdtLJK0PXEj1/HCl9iqllIgkAQcALVuk1jVtNboBcCBVKW7d6vpgu21xcxz7A7uV12cCC6mKrhERMYn0ukQo6VyqfDBD0gjwYaohPLF9KrAA2AdYTNVS9LBy6uuBlwGbS3pr2fZW29dT1V5uQVV7eT1VjWZbTatGLwIeAa4FljY8B6p/E3y3jEjzBdvzgC1Hi64lsz91AteLiIg+6fXIMrYP7rDfwBEttp9NNd51q3P2mGgcTRPhNrb3mujFgV3LuKRPBS6R9NOmJ5bWRYcDTNt009W4dURErK5hGnS7afeJH0p67kQvXvp3YPsBqvrcnYFf1upwZ1L1/2h17jzbc23PnbZRhjaNiOgrg1eoq2WqaJoIX0rVc/+20jdjkaQb250g6cmjM1SUQbr3pHpoOR84tBx2KFW1a0RETDaDbyzTF02rRvdejWtvCVxYGtasC3zV9nckXU01iPfbgV8Af7sa146IiJ5qNjrM2qBRIrT9c4DyrG/DhufcCTy/xfZfAy+fQIwRETEIU6hU142m3Sf2Az4JbEX1TG874Fbg2b0LLSIiBmZydKjvi6bPCD8C7AL8zPZsqhLdD3oWVUREDN6QPCNsmgj/XKo015G0ju3LgI5zPEVExFSmLpepoWljmYclbQRcQdVr/wFgWe/CioiIgZtCpbpuNC0R7k81vM17gO9QTXexb6+CioiISSBVoys5zvYK28tsn2n7ZDI+aETE2isd6lfxyhbbVqdvYURETBVDUiJs+4xQ0juBfwD+YsxIMtNJq9GIiLXbkHSf6NRY5qvAt4GTWHnewEdtP9SzqCIiYuA0hUp13WhbNWr7Edt3Ax8C7i8jzMwG3ixpkz7EFxERg9BttegUSqJNnxF+HVgu6RnAl6mS4Vd7FlVERAyYqqrRbpYpomkiXGF7GfBa4DO23wPM7F1YERExcENSImzaof7Pkg4GDgFeXbat15uQIiJiUphCyawbTUuEhwEvAU60fZek2cDZvQsrIiIGbkhKhI0Soe1bbB9l+9yyfpftj/Y2tIiIGBiDVqirpRNJp0l6QNJN4+yXpJMlLS6Twu9U23eopNvLcmht+wvL5PGLy7kdA2mbCCWdX34uKkGstHR8lxERMXX1vkR4BrBXm/17AzuU5XDgFABJmwEfBl4M7Ax8WNKm5ZxTyrGj57W7PtD5GeHR5WfGFY2IiDXK9hWSZrU5ZH/gLNsGrpS0iaSZwG7AJaP92SVdAuwlaSGwse0fle1nAQdQ9YcfV9tEaHtJ+fnzBu8pIiLWIpOgQ/3WwD219ZGyrd32kRbb2+o0xNqjtCng2t640w0iImKK6r4v4AxJ19TW59meN4HzWwXg1djeVqcS4XQASScA9wNfKTd6E9V4oxERsTZaMy0/H7Q9t4vzR4Bta+vbAPeV7buN2b6wbN+mxfFtNe0+8Srb/277Udu/tX0KcGDDcyMiYioafPeJ+cAhpfXoLsAj5ZHdxcCekjYtjWT2BC4u+x6VtEtpLXoIcFGnmzTtUL9c0puA86je3sHA8om/p4iImCp6/YxQ0rlUJbsZkkaoWoKuB2D7VGABsA+wmGpy+MPKvockfQS4ulzqhNpEEO+kao36JKpGMm0bykDzRPhG4LNlMdUUTG9seG5ERExFPU6Etg/usN/AEePsOw04rcX2a4DnTCSORomwzECx/0QuHBERU5cMWjHoKPqjaYkwIiKGzRSaQaIbSYQREdHa4PsR9kUSYUREtDQJOtT3RacO9f/Ybr/tT63ZcCIiYtJIIgTSaT4iYjg5JUIAbP9zvwKJiIhJJonwCZI2BN4OPBvYcHS77bf1KK6IiBi0IUmETYdY+wrwNOBVwOVU47c92qugIiJi8OTulqmiaSJ8hu3/C/ze9pnA/wKe27uwIiJi4AY/1mhfNO0+8efy82FJz6GaiWJWTyKKiIjBm2Klum40TYTzygjf/5dqNPCNgON6FlVERAxeEuETbH+pvLwc2L534URExKSRRPgESRtQzT84q36O7RN6E1ZERAySGJ6q0aaNZS6imn1iGfD72tKRpGmSrpP0zbJ+hqS7JF1fljmrE3hERPRYGsusZBvbe63mPY4GbgU2rm17r+0LVvN6ERHRa0PUWKZpifCHkibcXULSNlRdLb7U6diIiJhkhqRE2DQRvhS4VtJtkm6UtEjSjQ3O+wzwPmDs9I4nlut8ujx/jIiIyWZIEmHTqtG9J3phSfsCD9i+VtJutV3HUvVDXB+YB7wfWKXRjaTDgcMBpm266URvHxERXRqWGeoblQht/xzYBHh1WTYp29rZFdhP0t3AecAeks62vcSVpcDpwM7j3HOe7bm2507b6MkN305ERKwR3ZYGp1CJsFEilHQ0cA7w1LKcLeld7c6xfaztbWzPAg4Cvmf7zZJmlmsKOAC4qYv4IyKiR/ox1qikvcpjt8WSPtBi/3aSLi2P0xaWtidI2r3W++B6SX+UdEDZN6HeCU2rRt8OvNj278tNPgb8CPi3hufXnSNpC6puKtcD71iNa0RERK/1uFQnaRrweeCVwAhwtaT5tm+pHfYJ4CzbZ0raAzgJeIvty4A55TqbAYuB79bOa9w7oWkiFLC8tr68bGvE9kJgYXm9R9PzIiJicPrQfWJnYLHtOwEknUfVZ72eCHcE3lNeXwZ8o8V1Xgd82/ZjqxNE01ajpwNXSTpe0vHAlcCXV+eGERExRXT/jHCGpGtqy+Fj7rA1cE9tfaRsq7uBamQzgNcA0yVtPuaYg4Bzx2xr3Duh6Vijn5K0kKobhYDDbF/X5NyIiJiC1kyDlwdtz22zv1XN4ti7HgN8TtJbgSuAe6lGOasuULU7eS5wce2cRr0TRrVNhJI2tv3bUv96d1lG921m+6F250dExNQkJvD8a/WNANvW1rcB7qsfYPs+4LUAkjYCDrT9SO2Q1wMX2v5z7Zwl5eVSSadTJdNxdSoRfhXYF7iWlbO0ynpmooiIWFv1/hnh1cAOkmZTlfQOAt5YP0DSDOAh2yuoSnqnjbnGwWV7/ZyZtpc07Z3QNhHa3rf8nN3x7URExFql1x3qbS+TdCRVteY04DTbN0s6AbjG9nxgN+AkSaaqGj3i8fikWVQlysvHXHpCvRM6VY3u1OFN/KTd/oiImML60Cne9gJgwZhtx9VeXwC07AZh+25WbVwz4d4JnapGP1l+bgjMpWq9I+B5wFVUjWciImJtk9knKrZ3t7078HNgpzLk2QuBF1B1XoyIiLXVkAyx1rRD/bNsLxpdsX1TJtSNiFi7DUuJsGkivFXSl4CzqfL8m6km242IiLVVEuFKDgPeSTXbPFQtd07pSUQRETEppERYY/uPwKfLEhERa7sp9pyvG526T5xv+/WSFtHiI7H9vJ5FFhERg5VECDxRFbpvrwOJiIjJQwzPDPWdRpZZUuaL+rLtV/QppoiImATk4SgSdpyGyfZy4DFJT+lDPBERMRl024dwCuXQpq1G/wgsknQJ8PvRjbaP6klUERExcGk1urJvlSUiIoZFEuETbJ8p6UnA023f1uOYIiJiEhiWEmHHZ4QAkl5NNZXFd8r6HEnzexlYREQM2JA8I2yUCIHjgZ2BhwFsXw9kjsKIiLVVmX2im2WqaPqMcJntR6rJfh83hd5mRERM2JD8lW+aCG+S9EZgmqQdgKOAH/YurIiIGKSqQ/1wZMKmVaPvAp4NLAXOBX4LvLtXQUVExOClarTG9mPAB8sSERFruynW4KUbnQbdbtsy1PZ+azaciIiYLDLWaOUlwD1U1aFXUVUbR0TEMOhDiVDSXsBngWnAl2x/dMz+7YDTgC2Ah4A32x4p+5YDi8qhvxgtnEmaDZwHbAb8BHiL7T+NF0OnZ4RPA/4JeE4J9JXAg7Yvt335BN5rRERMMb1+Rlgmdfg8sDewI3CwpB3HHPYJ4Kwy7d8JwEm1fX+wPacs9RrKjwGftr0D8Bvg7e3iaJsIbS+3/R3bhwK7AIuBhZLe1fktRkTElGXA7m7pbGdgse07S4ntPGD/McfsCFxaXl/WYv9KVPXz2wO4oGw6Ezig3TkdW41K2kDSa4GzgSOAk4H/7HReRERMbWugRDhD0jW15fAxt9ia6vHbqJGyre4G4MDy+jXAdEmbl/UNy3WvlDSa7DYHHra9rM01V9KpscyZVNWi3wb+2fZN7Y6PiIi1SPfPCB+0PbfN/lbtTsbe9Rjgc5LeClwB3AuMJrmn275P0vbA9yQtoure1+maK+nUWOYtVNMu/SVwVG1kGQG2vXGH8yMiYgqS3Y8O9SPAtrX1bYD76gfYvg94LYCkjYADbT9S24ftOyUtBF4AfB3YRNK6pVS4yjXH6vSMcB3b08uycW2ZniQYEbF260OH+quBHSTNlrQ+cBCwUrc9STMkjeaqY6lakCJpU0kbjB4D7ArcYttUzxJfV845FLioXRBNR5aJiIhhs7qzTjTsjF9KbEcCFwO3AufbvlnSCZJGW4HuBtwm6WfAlsCJZftfAddIuoEq8X3U9i1l3/uBf5S0mOqZ4ZfbxdF0rNGIiBgy/RgmzfYCYMGYbcfVXl/AEy1A68f8EHjuONe8k6pFaiNJhBERsSoDQzLodhJhRES0Nhx5MIkwIiJam0ozSHSj541lJE2TdJ2kb5b12ZKuknS7pK+VlkIRETHZ9H5kmUmhH61Gj6ZqDTRqQmPARUTEYAzLfIQ9TYSStgH+F/Clsj7hMeAiIqL/5GqG+m6WqaLXzwg/A7wPmF7WG48BV8akOxxg2qab9jjMiIhYxZDMR9izEqGkfYEHbF9b39zi0Jb/bLA9z/Zc23OnbfTknsQYERHjk93VMlX0skS4K7CfpH2ADYGNqUqIExoDLiIiBqDh6DBrg56VCG0fa3sb27Ooxo/7nu03McEx4CIiYhC6bDE6hUqEgxhrdEJjwEVExGAMS6vRvnSot70QWFheT2gMuIiIGJApVKrrRkaWiYiIVRk0JK1GkwgjIqK1lAgjImKYTaVO8d1IIoyIiNZSIoyIiKFlhmZkmSTCiIhYhZhao8N0I4kwIiJaSyKMiIihNiSJcBAjy0RExGQ3+oywm6UBSXtJuk3SYkkfaLF/O0mXSrpR0sIyvR+S5kj6kaSby7431M45Q9Jdkq4vy5x2MaREGBERLfX6GaGkacDngVdSTct3taT5tm+pHfYJ4CzbZ0raAzgJeAvwGHCI7dslbQVcK+li2w+X895r+wIaSIkwIiJa6/2g2zsDi23faftPwHnA/mOO2RG4tLy+bHS/7Z/Zvr28vg94ANhidd5mEmFERKzKhhUrultghqRrasvhY+6yNXBPbb3VZO03AAeW168BpkvavH6ApJ2B9YE7aptPLFWmn5a0Qbu3mqrRiIhorft+hA/anttmf5PJ2o/RwQ44AAANgUlEQVQBPifprcAVwL3AsscvIM0EvgIcans04mOB+6mS4zyqWY9OGC+IJMKIiGipD/0IR4Bta+urTNZeqj1fCyBpI+BA24+U9Y2BbwEfsn1l7Zwl5eVSSadTJdNxpWo0IiJa6/0zwquBHSTNlrQ+1STu8+sHSJohaTRXHQucVravD1xI1ZDmP8acM7P8FHAAcFO7IJIIIyJiVQZWuLul0y3sZcCRwMXArcD5tm+WdIKk/cphuwG3SfoZsCVwYtn+euBlwFtbdJM4R9IiYBEwA/iXdnGkajQiIlpoXKrr7i72AmDBmG3H1V5fAKzSDcL22cDZ41xzj4nEkEQYERGtDcnIMkmEERHRWhJhREQMrdFnhEMgiTAiIlowrFg+6CD6IokwIiJWlRJhREQMvTwjjIiIoZZEGBERw6s//QgngyTCiIhYlRmdQWKtl0QYERGtpUQYERFDLYkwIiKGV7OBs9cGSYQREbEqg5enQ31ERAyzVI1GRMTQstNqNCIihlxKhBERMcycEmFERAyvjCwTERHDLLNPRETE0PNwVI2u06sLS9pQ0o8l3SDpZkn/XLafIekuSdeXZU6vYoiIiNVjwCvc1dKEpL0k3SZpsaQPtNi/naRLJd0oaaGkbWr7DpV0e1kOrW1/oaRF5ZonS1K7GHpZIlwK7GH7d5LWA/5b0rfLvvfavqCH946IiG7YPe9QL2ka8HnglcAIcLWk+bZvqR32CeAs22dK2gM4CXiLpM2ADwNzqfL2teXc3wCnAIcDVwILgL2AbzOOnpUIXfldWV2vLMNR4RwRsTbwiu6WznYGFtu+0/afgPOA/cccsyNwaXl9WW3/q4BLbD9Ukt8lwF6SZgIb2/6RbQNnAQe0C6KnzwhLtr8WeAbwedtXSXoncKKk46je3AdsL21x7uFUGR1g6V1HH3NTL2Ndg2YADw46iIamUqwwteJNrL2RWHvnmfWVR/nNxf/lC2Z0ec0NJV1TW59ne15tfWvgntr6CPDiMde4ATgQ+CzwGmC6pM3HOXfrsoy02D6uniZC28uBOZI2AS6U9BzgWOB+YH1gHvB+4IQW584r+5F0je25vYx1TUmsvTOV4k2svZFYe2dMwsL2Xv24bYttY2sOjwE+J+mtwBXAvcCyNuc2ueZKelY1ulIE9sPAQmAv20tKtelS4HSqonFERAyfEWDb2vo2wH31A2zfZ/u1tl8AfLBse6TNuSPl9bjXHKuXrUa3KCVBJD0JeAXw01J/S2nFcwAwVao8IyJizboa2EHSbEnrAwcB8+sHSJohaTRXHQucVl5fDOwpaVNJmwJ7AhfbXgI8KmmXkmcOAS5qF0Qvq0ZnAmeW54TrAOfb/qak70nagqr4ej3wjgbXmtf5kEkjsfbOVIo3sfZGYu2dvsdre5mkI6mS2jTgNNs3SzoBuMb2fGA34CRJpqoaPaKc+5Ckj1AlU4ATbD9UXr8TOAN4ElVr0XFbjALIQzKETkRERCt9eUYYERExWSURRkTEUJs0iVDS35ah2FZImlvbPkvSH2pDsp06zvmbSbqkDLVzSXl42u9YXynp2jK0z7VlFIRW5x8v6d7ae9qn37GWfceWIYhuk/Sqcc6fLemq8rl+rTzQ7rlyr9HP525J149z3N3l875+bPPvfmr6O+00nFQ/SPpXST8tQ1ZdONqorcVxA/tsO31OkjYo35HF5fs5q5/x1eLYVtJlkm4t/58d3eKY3SQ9UvtuHDeIWEssbX+nqpxcPtcbJe00iDj7zvakWIC/ourQuRCYW9s+C7ipwfkfp+qcD/AB4GMDiPUFwFbl9XOAe8c5/3jgmAF/rjtSdVTdAJgN3AFMa3H++cBB5fWpwDsH8N34JHDcOPvuBmb0O6bV+Z1SNQa4A9ieqh/tDcCOA4h1T2Dd8vpj4/2/MqjPtsnnBPwDcGp5fRDwtQH93mcCO5XX04GftYh1N+Cbg4hvor9TYB+qhiUCdgGuGnTM/VgmTYnQ9q22b+viEvsDZ5bXZ9JhSJ1ujBer7etsj/ZXuZlqVIUNehVHE20+1/2B82wvtX0XsJgxfTpL0+M9gNFxYXv6ubZSYng9cG4/79sjTYaT6jnb37W9rKxeycp9riaDJp9T/f/3C4CXl+9KX7nqF/2T8vpR4FY6jGIyye1PNa6nbV8JbDLa5W1tNmkSYQezJV0n6XJJfzPOMVu66j9C+fnU/oXX0oHAdW4xfFxxZKl6OK2X1bhtjDc8Ud3mwMO1P5odhyrqgb8Bfmn79nH2G/huqYo+fJxj+qXT77TJZ95vb2P8puWD+mybfE6PH1O+n49QfV8HplTPvgC4qsXul6iaiefbkp7d18BW1ul3Ohm/oz3X1/kIJf0X8LQWuz5oe7wOj0uAp9v+taQXAt+Q9Gzbv+1ZoKx2rKPnPpuqymnPcQ45BfgI1ZfyI1RVf2/rc6xNhiGa8FBFE9Ew7oNpXxrc1fZ9kp4KXCLpp7avWFMx1rWLl2a/055+nivdqMFnK+mDVENVnTPOZfr22Y4x8O/mREnaCPg68O4Wf5t+AmznaiaefYBvADv0O8ai0+90Un2u/dLXRGj7FatxzlKqKZ2wfa2kO4C/BMY+6P2lpJm2l5Si/AP9jhVA1VxZFwKH2L5jnGv/snb8F4FvrlaQT1xvdWLtOLQR1YDBm0hat/yru+NQRRPRKW5J6wKvBV7Y5hr3lZ8PSLqQqlqtJ3+sm37ObX6nTT7zNaLBZ3sosC/wcpeHQy2u0bfPdowmn9PoMSPle/IU4CEGQNU0c18HzrH9n2P31xOj7QWS/l3SDNt9H5C7we+0b9/RyWTSV42qGqptWnm9PdW/pO5sceh8YHRixkPpMKROL5TWd98CjrX9gzbH1evcX8NghpmbDxxUWt/Npvpcf1w/oPyBvAx4XdnU78/1FcBPbY+02inpyZKmj76mKoEPZMi+hr/TjsNJ9YOkvagGu9/P9mPjHDPIz7bJ51T///11wPfGS+i9VJ5Lfhm41fanxjnmaaPPLyXtTPV399f9i/LxOJr8TucDh5TWo7sAj4w+clqrDbq1zuhC9cdjhKr090uqMeOgetZ2M1XLsZ8Ar66d8yVKS0iq5wOXAreXn5sNINYPAb+nGjpudHlqi1i/AiwCbqT64s3sd6xl3wepWufdBuxd276AJ1q/bk+VIBcD/wFs0MfvxBnAO8Zs2wpYUIvthrLcTFXtN6jvb8vfaT3esr4PVcvCOwYVb/ld3lP7jp46NtZBf7atPieqWWr2K683LN/HxeX7uf2APsuXUlUd3lj7PPehGjryHeWYI2t/w64E/npAsbb8nY6JVVQT5d5Rvs9zBxFrv5cMsRYREUNt0leNRkRE9FISYUREDLUkwoiIGGpJhBERMdSSCCMiYqglEUZExFBLIoxJSdJratPWjC4rJO29hu9zhqTXdThmlqSbyuu5kk7ucOwb12SMTak2ZdkEz3tDmXanqxGOIqaqJMKYlGxfaHvO6AL8O/B94OIBx3WN7aPaHDILGEgiLO4on1djtr8G/O8exRMx6SURxqQn6S+B44C32F5Rhn/6V0k3qZpk9A3luN0kLZR0gaqJZ8+pDW11nKSryznzRre3uecLy2wBPwKOqG3fbbTkJOl/1kqr15Xhqz4K/E3Z9p5SSvu+pJ+U5a8bxPoiST8s9/+xpOmSppX3fLWqGS7+vsHnNqtc+0vlfZ8j6RWSfqBqouWdO10jYhgkEcakVgY0/irVpLe/KJtfC8wBnk81Hum/1sb6fAHwbqqJh7cHdi3bP2f7RbafAzyJasDpdk4HjrL9kjbHHAMcUUpgfwP8gWpS6O+XkuynqQZ/f6XtnYA3APVq1VViLWNrfg042vbo+/sD8HaqcR9fBLwI+LsyRmwnzwA+CzwPeBZVafWlJfZ/anB+xFoviTAmu48AN9s+r7btpcC5tpe7msnjcqrkAPBj2yO2V1CN+zirbN9d0lWSFlFNNjzunHCSngJsYvvysukr4xz6A+BTko4qxy9rccx6wBfLff+DKumNahXrM4Eltq+GauaCct09qQZDvp5qvrvNaTaVz122F5V73Axc6mpcxUU88dlEDLW+TsMUMRGSdqMadH2nsbvanFafCHk5sK6kDameMc61fY+k46kGbR731jSYg832RyV9i2qQ5SsltZr66D1Ug50/n+ofnn9sF2ubewt4l+2JPiOt32NFbX0F+f8/AkiJMCYpVTO8n041r+OjY3ZfAbyhPDfbAngZY6aQGmM06T2oagLVtq1EbT8MPCLppWXTm8aJ8S9KaetjVPNjPgt4FJheO+wpVCW8FcBbgGnt7g38FNhK0ovKPaarmm/vYuCdpaoYSX9ZptKJiC7lX4QxWb0DeCpwyph2LScB5wMvoZpOxsD7bN8v6VmtLmT7YVWT5S4C7qaa766Tw4DTJD3G+C1V3y1pd6rS3C3At6lKWssk3UA1jdS/A1+X9LdUczv+vt1Nbf+pNP75N0lPono++AqqabxmAT8pjWp+BRzQ4H1ERAeZhiliLSFpFvDN0iBooufuRtUgqVMjooi1TqpGI9Yey4GnrE6HeqqS6296ElXEJJcSYUREDLWUCCMiYqglEUZExFBLIoyIiKGWRBgREUPt/wNUuPOqGPTR2wAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "U: 1.000000 \n", "V: 1.000000 \n", "temp: 20.000000 \n" ] } ], "source": [ "fieldset_flat = FieldSet.from_data(data, dims, mesh='flat')\n", "fieldset_flat.U.show()\n", "for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temp]:\n", " print('%s: %f %s' % (fld.name, fld[0, 0, 40, -5], fld.units))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Indeed, in this case all Fields have the same default UnitConverter object. Note that the coastlines have also gone in the plot, as .show() recognises that this is a flat mesh." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### UnitConverters for Diffusion fields" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The units for Brownian diffusion are in $m^2/s$. If (and only if!) the diffusion fields are called kh_zonal and kh_meridional, Parcels will automatically assign the correct Unitconverter objects to these fields." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Kh_zonal: 1.380091e-08 \n", "Kh_meridional: 8.098704e-09 \n" ] } ], "source": [ "kh_zonal = 100 # in m^2/s\n", "kh_meridional = 100 # in m^2/s\n", "\n", "fieldset.add_field(Field('Kh_zonal', kh_zonal*np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid))\n", "fieldset.add_field(Field('Kh_meridional', kh_meridional*np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid))\n", "\n", "for fld in [fieldset.Kh_zonal, fieldset.Kh_meridional]:\n", " print('%s: %e %s' % (fld.name, fld[0, 0, 40, -5], fld.units))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, the unitconverters are GeographicPolarSquare and GeographicSquare, respectively.\n", "\n", "Indeed, multiplying with $(1852\\cdot60)^2$ returns the original value" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100.0\n" ] } ], "source": [ "deg_to_m = 1852*60\n", "print(fieldset.Kh_meridional[0, 0, 40, -5]*deg_to_m**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding a UnitConverter object to a Field" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, to summarise, here is a table with all the conversions\n", "\n", "| Field name | Converter object | Conversion for mesh='spherical'| Conversion for mesh='flat'|\n", "|-------|-----------------|-----------------------------------|\n", "| 'U' | GeographicPolar| $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$ | 1 |\n", "| 'V' | Geographic | $1852 \\cdot 60$ | 1 |\n", "| 'Kh_zonal' | GeographicPolarSquare | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1 |\n", "| 'Kh_meridional' | GeographicSquare | $(1852 \\cdot 60)^2$ | 1 |\n", "| All other fields | UnitConverter | 1 | 1 |\n", "\n", "Only four Field names are recognised and assigned an automatic UnitConverter object. This means that things might go very wrong when e.g. a velocity field is not called U or V. \n", "\n", "Fortunately, you can always add a UnitConverter later, as explained below:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0\n" ] } ], "source": [ "fieldset.add_field(Field('Ustokes', np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid))\n", "print(fieldset.Ustokes[0, 0, 40, -5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This value for Ustokes of course is not as expected, since the mesh is spherical and hence this would mean 1 degree/s velocity. Assigning the correct GeographicPolar Unitconverter gives" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.1747725785927632e-05\n", "0.9999999999999999\n" ] } ], "source": [ "from parcels.tools.converters import GeographicPolar \n", "\n", "fieldset.Ustokes.units = GeographicPolar()\n", "print(fieldset.Ustokes[0, 0, 40, -5])\n", "print(fieldset.Ustokes[0, 0, 40, -5]*1852*60*np.cos(40*np.pi/180))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, the UnitConverter can be set when the FieldSet or Field is created by using the fieldtype argument (use a dictionary in the case of FieldSet construction." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.1747725785927632e-05\n" ] } ], "source": [ "fieldset.add_field(Field('Ustokes2', np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid, fieldtype='U'))\n", "print(fieldset.Ustokes2[0, 0, 40, -5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using velocities in units other than m/s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some OGCM store velocity data in units of e.g. cm/s. For these cases, Field objects have a method set_scaling_factor(). \n", "\n", "If your data is in cm/s and if you want to use the built-in Advection kernels, you will therefore have to use fieldset.U.set_scaling_factor(100) and fieldset.V.set_scaling_factor(100)." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0\n" ] } ], "source": [ "fieldset.add_field(Field('Ucm', 0.01*np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid))\n", "fieldset.Ucm.set_scaling_factor(100)\n", "print(fieldset.Ucm[0, 0, 40, -5])" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }