{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Computing water masses volume with f2py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tutorial for LDEO : pick if you are running on byrd (True) or not (False)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "byrd=True" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# get the World Ocean Atlas annual T/S data\n", "\n", "if not byrd:\n", " # wget in disguise in the Makefile\n", " !make getdata\n", " woadir='./'\n", "else:\n", " # if running on byrd, we are gonna use files stored in ~rdussin/demo_files/\n", " woadir='~rdussin/demo_files/'" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# we need to compile the fortran code into something python can use\n", "! source deactivate pangeo; make > /dev/null; source activate pangeo" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import lib_watermass\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "#---------------------------------------------------------------\n", "# sample data set (use make getdata to download WOA annual data)\n", "fnameT = woadir + 'woa13_5564_t00_01.nc'\n", "fnameS = woadir + 'woa13_5564_s00_01.nc'\n", "\n", "woa = xr.open_mfdataset([fnameT,fnameS],decode_times=False)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "dask.array\n", "Coordinates:\n", " * lat (lat) float32 -89.5 -88.5 -87.5 -86.5 -85.5 -84.5 -83.5 -82.5 ...\n", " * lon (lon) float32 -179.5 -178.5 -177.5 -176.5 -175.5 -174.5 -173.5 ...\n", " * depth (depth) float32 0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0 ...\n", " * time (time) float32 6.0\n", "Attributes:\n", " standard_name: sea_water_temperature\n", " long_name: Objectively analyzed mean fields for sea_water_temperatur...\n", " cell_methods: area: mean depth: mean time: mean\n", " grid_mapping: crs\n", " units: degrees_celsius" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "woa.t_an" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "#---------------------------------------------------------------\n", "# compute metrics (dx,dy,dz) of ocean grid cells and mask\n", "Rearth = 6378e+3\n", "nx = woa.lon.shape[0]\n", "ny = woa.lat.shape[0]\n", "nz = woa.depth.shape[0]\n", "\n", "dxflat = (2*np.pi*Rearth/360) * (woa.lon_bnds[:,1] - woa.lon_bnds[:,0])\n", "dyflat = (2*np.pi*Rearth/360) * (woa.lat_bnds[:,1] - woa.lat_bnds[:,0])\n", "dzcolumn = (woa.depth_bnds[:,1] - woa.depth_bnds[:,0]).values\n", "\n", "dxflat2d, dyflat2d = np.meshgrid(dxflat,dyflat)\n", "lon2d, lat2d = np.meshgrid(woa.lon,woa.lat)\n", "\n", "dx = dxflat2d * np.cos(2*np.pi*lat2d/360)\n", "dy = dyflat2d\n", "\n", "dz = np.empty((nz,ny,nx))\n", "for jj in np.arange(ny):\n", " for ji in np.arange(nx):\n", " dz[:,jj,ji] = dzcolumn" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "#---------------------------------------------------------------\n", "# compute mask from missing value\n", "mask = np.ma.array(np.ones(woa.t_an.squeeze().shape))\n", "mask[np.isnan(woa.t_an.values.squeeze())] = 0" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home2/nhn2/miniconda3/envs/pangeo/lib/python3.6/site-packages/IPython/core/magics/pylab.py:160: UserWarning: pylab import has clobbered these variables: ['plt']\n", "`%matplotlib` prevents importing * from pylab and numpy\n", " \"\\n`%matplotlib` prevents importing * from pylab and numpy\"\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAApkAAAKvCAYAAAAhjNTiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X3QbVldH/jvLxcB8Q2hwek37DZpR1vGq+SmZUKNQSHpxnFop0qtZkZtCZOeOOBb4gvEKTGOVKEmYciM0bmRDpAyNogmdBniVRGKmanw0iiNAqN0gMD1EhsQjFOMkL53zR/PuXL69DnPedv7Ofuc/flUPXWfs89+9lln7bV3r16/316rWmsBAIAu/YVdFwAAgMOjkwkAQOd0MgEA6JxOJgAAndPJBACgczqZAAB0bmkns6rurKr7q+r3prZ9VVW9uareUVX3VNVNk+1VVf+4qu6rqndW1ZP6LDwAANuZ19ebeX+j/t0qI5kvT3LLzLafSvL3W2tfleRHJ6+T5BlJbpj83JHkZ1cpBAAAO/PyPLSvN22j/t3STmZr7U1J/nh2c5LPn/z+BUkuTH6/Nckr25E3J3l0VV25SkEAADh5C/p60zbq3z1sw/J8X5JzVfUPctRR/auT7Vcn+dDUfucn2z48e4CquiNHveF8zqPqL3/ZX3r4hkUBAHiot7/zUx9trT1ul2W4+es+p33sjy/usgh5+zs/9a4kfza16Wxr7ewah1i5fzdt007mdyX5/tbaL1fVtyZ5WZKnJ6k5+85dt3Ly5c4myZnTj2xvPfeEDYsCAPBQp65877/fdRk+9scXs+s+zqkr3/tnrbUzWxxi5f7dtE2fLr89ya9Mfv+lJDdNfj+f5Nqp/a7JZ0LpAADsn436d5t2Mi8k+WuT378+yXsnv9+d5DsmTyE9OcmftNaOHUoFAGDQNurfLQ2XV9UvJnlqkiuq6nySFyb5W0leWlUPy1GM/47J7q9L8g1J7kvyySTP3uCLAAAchJbkUi7tuhjHWtDX+6wkaa39XDbs3y3tZLbWnrXgrb88Z9+W5LmrfDAAALt3TF/v8vsb9e82ffAHAIClWi62YY9k9sWykgAAdE4nEwCAzgmXAwD05OjBn6VTSh4kI5kAAHTOSCYAQI+GPoVRX4xkAgDQOZ1MAAA6J1wOANCTlpaLzYM/AADQCSOZAAA9MoURAAB0RCcTAIDOCZcDAPSkJbkoXA4AAN0wkgkA0CMP/gAAQEeMZMLEzVedzrkL9+66GA9y81Wnd12EwdUJAPtBJxMAoCctseIPAAB0xUgme6fPEPKmx54OKa9yjHMX7h1EKHwVi8opjA6wmku7LsCOGMkEAKBzOpkAAHROuBwAoCctbbQr/uhksnP7kpt4nHW/wxi/cyKPc2hWzbe9vN/Qz9+iaciGOD0ZjIFwOQAAnTOSCQDQl5ZcHGe03EgmAADdM5LJThxCTiKr2ddcuC7mTD3umF3VSx/zsi7ad+i5jcd9zz7qHlbRYp5MAADojE4mAACdEy5nK6uEoPZl+pNF07awmqGf38t2eV77qKNVj7kv52cbq57befvNhtrHUF+clMrF1K4LsRNGMgEA6JyRTACAnrQkl0xhBAAA3TCSyYOsuszc7Pabrzq99G9n31+0/STM+2w5mJsbY/7aGL/z0K16TuZd67Pb5uVnDn0KJxganUwAgB558AcAADpiJHNE1g2FL3t/0XQf0yGl2X0Whadn/37Rtk0dF+4XJt/ePoUR152qalE7NN3N/lrlnrbO38FxWoxkAgBAZ3QyAQDonHA5AECPLrVxhsurtd3PEHrm9CPbW889YdfFOChd5jAel3u2L/mMq+RS7ct3GaJDyVWTZ/kZ614PQ6yvVZaP7NoQ62HMTl353re31s7ssgxf8ZUPb//iV79ol0XIV33x+Z3Ug5FMAICeePAHAAA6pJMJAEDn5GSOyJBzDvtYYnJXuVFDruc+yUXbP9u01b7P97K5V4d8nQ3hWtikfla5Dw/hu61jCDmZX/6Vj2iv/NUrd1mE3PTF/34n9WAkEwCAznnwBwCgR2Odwkgn88CsGiKZF/LoO0SyaLnJPuwypNPlcphDt2+hM46s0i53fW43XfqRI+veh6b3H9M9jH4JlwMA0DkjmQAAPTFPJgAAdMhI5kgtWm5tk7+b97eL9ltn+zbLv236/Vid+txPq17DXRz/uLzKQ8653DTnddnfbXqOtsmxdJ2zDZ1MAIDeVC62cQaOx/mtAQDolZFMAICetCSXRjqmp5N5QLbNZTrJv1+Ub9nn/Gyr5Ip1bZu80qGSo7V/+s7DXPQZs7mXl/c5xOtiXZt8/+n66/rz9mHuVPbPOLvWAAD0ykgmAECPxjpP5iA6mX/wzkftJJQ567hpNfbdSYSn1pmSY3b7Sde7KY5g+3veKveURWF0unGSdbnt9FQn6TNlfe9OyzF2g+hkAgAcotZMYQQAAJ3RyQQAoHPVWjt+h6o7k3xjkvtba0+c2v7dSZ6X5IEk/7q19kOT7S9I8pwkF5N8T2vt3LJCfH49pn3iw1ds/CVWtYucxL70OZVFlzZZOq7PaYzW0fX5PJQ8tCHkW52k6Ryzbc7hUOrtuOmM1r2vrDPtzaG0fza3i/9G/mZ7zdtba2dO7IPn+NL/4rPb/3b39bssQm75kvfspB5Wycl8eZL/PckrL2+oqq9LcmuSr2ytfaqqHj/ZfmOS25J8RZKrkvxmVX1pa+1i1wUHAGC4lnYyW2tvqqrrZjZ/V5IXt9Y+Ndnn/sn2W5PcNdn+/qq6L8lNSf5tZyUGANgTLcnFkWYnbvp0+Zcm+a+q6kVJ/izJD7TW3pbk6iRvntrv/GTbQ1TVHUnuSJInXP2wJP2Gy08iVLPuZ3QRFt7k7056tY11P2s6ZLfOtEjzbBOq63tarX0KIw4l1HtSZs9JF+do1RSRrj5j0TFX+Syrv9CHfbjX0a1Nu9YPS/KFSZ6c5AeTvLqqKpk72+jcpM/W2tnW2pnW2pnHPfbUhsUAAGCINh3JPJ/kV9rRU0NvrapLORqKPJ/k2qn9rklyYbsiAgDsK/NkrutfJfn6JKmqL03y8CQfTXJ3ktuq6hFVdX2SG5K8tYuCAgCwP5aOZFbVLyZ5apIrqup8khcmuTPJnVX1e0k+neT2yajmu6rq1UnenaOpjZ67ypPlfSwruQ+5H7sq40l8bh/ncJVjDmX6o3lmyza08k0bU77dEM7DcdMKrbL/5f2m8637/l6rHn8I9Qu71JJc8uDPfK21Zy1469sW7P+iJC/aplAAAOy3cXatAQDo1aYP/gAAsIKLbd7kO4dvkJ3MdfOTFu0rF2h3+qj7beYh7euz1/2MIbTJMeVb7rt1c9WH0L5m7dN8sEC3hMsBAOjcIEcyAQAOQUtZVnIfzAu3rDrFB4fvuHDi7HtdtY/pJTCHbNvyrRK27XsZzj7s033iuLIO+XsMuWxAv/aqkwkAsG8uWfEHAAC6oZMJAEDnDipcLvdnnHY5tcu6uYcnOZ1LV3mRs2VdpexDz88c0r3iuPoZUjmBzbRktA/+jPNbAwDQq4MayQQAGJKWGu2KP0YyAQDo3CBGMr/0Kz+Zc+dWn+9yUQ7T9Ha5TONxXP5fH+2gixzDLpc97TvncZuyDikfs4/67uKY7lXAoRpEJxMA4FBdGmngeJzfGgCAXg1+JHNeuG1ReGl6X6Hzw7DueTxumcd12hIPts7US4cYIu/7mMDhai25aMUfAADohk4mAACdG3y4HABgf1UuZZzzZO5lJ3PVnC+5U/tplSmqLls0fdHNV51eeTqjk1zqcd7nbvPZJ53/2NfnHZdLu8lxuiS/G2Aze9nJBADYBy0e/AEAgM4c7EimsNZ+2TRM2kX4dp1QbV/tat2Q/ZCmCerCEL/PEMsEPNiy6/TUlSdUEOY62E4mAMAQXBxp4Hic3xoAgF4ZyQQA6ElL5VIzhdHgrZI7Jxdz//SR+7bJMRdNh3SSzl24VxueYxfnZvp+45zAahZdn9PTyu3LlG1sT7gcAIDO7dVIJgDAvvHgDwAAdGSvRjJXycdYtPTgNnlVi/52evsu51ikf9PnrqslEBdZZ/nMsdj1d152ToDlFv13lMO1V51MAIB90pJcsqwkAAB0o1pruy5Dzpx+ZHvruSf0/jnbhLf6GNIXbnuoLuv5JOpXqKcby6YLWqeeF533rpbudN1yKE7i/rXK9dhnOU5d+d63t9bO9PYBK7j2iV/QvveXnrzLIuQHb/z1ndSDkUwAADqnkwkAQOc8+AMA0JMxP/gzqk7m0HKspssjz+tIl1P1nET99j2d0dgsOk/r5G0tmm5onWUiZ/dxfjkkfbfnVe637p3jMKpOJgDASbuY2nURdmKc47cAAPRKJxMAgM6NNly+bO6u6X1OImdk3Tn8xmC2LjaZK/Gk8n36nuvt0C1r/5vWbddz47o+2VcneS/kwVqr0T74M85vDQBAr0Y7kgkAcBIujnQkczSdzE3CcLsIf64yvcpYrTLtzOzvJ1mXpuTYzKrTnXRVr12H5GHodrl85DzHTSfoOjws4+xaAwDQq9GMZAIAnLSW5JJ5MgEAoBsHPZK5zRQ4uzSdT2hqo83ssr4sR7jcJuenr5xXy7ty6E7intTFf6sO935Ze/HgT1XdkuSlSU4l+fnW2otn3n9CklckefRkn+e31l533DGH/60BAOhNVZ1K8jNJnpHkxiTPqqobZ3b7n5O8urX21UluS/JPlh1XJxMAYNxuSnJfa+19rbVPJ7krya0z+7Qknz/5/QuSXFh20IMKl8+GBA5h6F34bnXT6QVDmQpK6Hw9q4Suj5vmZJswPIzBLqfwGus9sCW51Hb+4M8VVXXP1OuzrbWzU6+vTvKhqdfnk3zNzDF+LMmvV9V3J/mcJE9f9qEH1ckEAOAhPtpaO3PM+/N6wW3m9bOSvLy19g+r6r9M8s+r6omttUuLDqqTCQDQo4vDz048n+TaqdfX5KHh8OckuSVJWmv/tqoemeSKJPcvOujgvzUAAL16W5Ibqur6qnp4jh7suXtmnw8meVqSVNWXJ3lkko8cd9DBj2Suk9N2aPke03mFpllhLFZdZjI5vGseDk0f16jlJ7vXWnugqp6X5FyOpie6s7X2rqr68ST3tNbuTvJ3k/zTqvr+HIXSv7O1NhtSf5DBdzIBAPZVSw3hwZ+lJnNevm5m249O/f7uJE9Z55jC5QAAdE4nEwCAzg0yXD6dbzX2vIt5338oc0AOkXrZP9te4845DMeiZwm6/oxlPnNfeG8vZVjXpZGO6Y3zWwMA0KtBjmQCAByC1pKLe/DgTx+WdjKr6s4k35jk/tbaE2fe+4EkP53kca21j1ZVJXlpkm9I8skcPd7+2+sWauwhcvabdIbltrnGp+tXPcP6Ll9/XV4/Q1tC93IZTl2544KM3Crh8pdnMsP7tKq6Nslfz9HknJc9I8kNk587kvzs9kUEAGDfLB3JbK29qaqum/PWS5L8UJLXTm27NckrJ5NzvrmqHl1VV7bWPtxFYQEA9s0+zJPZh40e/KmqZyb5w9ba7Jj41Uk+NPX6/GTbvGPcUVX3VNU9H/nYxU2KAQDAQK394E9VPSrJjyT5G/PenrNt7pJDrbWzSc4myZnTjzx2WSIeqo+cGja33pQaw8hZGqp5dTNdd9o8dG/bpRrd0xY7WvFnnJP5bPJ0+V9Mcn2Se4+e88k1SX67qm7K0cjltVP7XpPkwraFBABgv6zdtW6t/W5r7fGttetaa9flqGP5pNbaf0hyd5LvqCNPTvIn8jEBAMZnlSmMfjHJU5NcUVXnk7ywtfayBbu/LkfTF92XoymMnt1ROQEA9tLFudmEh2+Vp8ufteT966Z+b0meu32xWNV0HoxctZOxLPfIediO+uuG+VpZhVxK+mTFHwCAnrSYwggAADpjJPOACI/tj22nCzlkpuda36K2tGw6KMbHfYeTpJMJANCb8c6TOc5vDQBAr4xkAgD06JIpjIBVzFsecjrH8ricp8t/Ky9quTFOz3USedVyt8dp0T3HPYk+CZcDANA5I5kAAD1pLbk40nkydTJhC+uGHYWkNnMIId5Vwv+bfEdtimVm28i8lJ9FTLfGNoTLAQDonJFMAIAemScTAAA6YiQTtiBXab+d1JRBx33ucflyq1iUM7fvOaxsb9lyozdfdXppO3SP215L5dJIH/wxkgkAQOd0MgEA6JxwOQBAjywrCfw5+Un9O4mcwWVzU87mM3ZdpuPay7J8uW3Lcwhzi7K5ZctIHrf93IV73evohE4mAEBPWuLBHwAA6IqRzAMiNLa9VUKYDM86S+PNXifTqRCL9um6POseZ1l5Fi0T2NX3YbiOmw5rk2NBl3QyAQB6ZMUfAADoiJFMAIC+tPGu+KOTeUBMWcIYrZtHtug6WZTXuMk1te3UQ13rc5omTt5sG9m0zcjBpG/C5QAAdM5IJgBAT1rGu+KPkUwAADpnJJPRk5e0f1Y9Z4vyD1dZbnJXtv3sRXMmysWE3Rnrgz9GMgEA6JxOJgAAnRMuZ5SEyPfXvDDwNudzUeh8lWl/hjg10LKlUYdSzrFb93z0cc/a9tphNS3C5QAA0BkjmQAAPTKSCQAAHTGSeWDkXS0nB2m+fc3POq7M63yfRfmY069nr6tF0wXBcabby6K2M6+tbdPO5v03QbulbzqZAAA9aSnhcgAA6IqRTCDJQ6fjOaRQ2qLvMxtCXPc7T4fRh5KqIoQ/fKtMiXXuwr1bp2fsui2CTiYAQI8uRbgcAAA6YSQTAKAvbbzzZOpkHqh5+TywqpPI51slR3LRvpu07Xn5bJvkYO7LVDCu//3R15KRJ/l5MI9wOQAAnTOSCQDQk5bxhsuNZAIA0DkjmQdMXib7ZtUl9rax7Rygs/OJzm5f5dizf+c6HZ8u8yK1n+EzkgkAAB3RyQQAoHPC5YyGaTuG6yTDfV22g67C7ozPJst/LkrVWPdvOVktJVwOAABdMZIJANCjZiQTAAC6YSTzwK2zdN8hkoc0fCfdHredwqgL2+TWMR7bTOm16zYOiU4mAECvLkW4HAAAOmEkEwCgJ62Nd8UfncwRGssydnKSTta68/7tsg0u+ux1c+C0Mba1To7wvHaoDTJkwuUAAHTOSCYAQI/GOk+mTiYHRehod9ap+6Gma3QV8jdFEetYtd1dfl+bYl/oZAIA9Mba5QAA0BmdTAAAOrc0XF5Vdyb5xiT3t9aeONn200n+mySfTvLvkjy7tfaJyXsvSPKcJBeTfE9r7VxPZWcL83J7xjK1ESyz7XXgOqJPcs/3z1gf/FllJPPlSW6Z2fYbSZ7YWvvKJH+Q5AVJUlU3JrktyVdM/uafVNWpzkoLAMBeWNrJbK29Kckfz2z79dbaA5OXb05yzeT3W5Pc1Vr7VGvt/UnuS3JTh+UFAGAPdPF0+d9M8qrJ71fnqNN52fnJtoeoqjuS3JEkT7jaQ+5DMOQQ32x4aHrKjyGXm/HZdZs0zQ0MS8t4l5Xc6sGfqvqRJA8k+YXLm+bs1ub9bWvtbGvtTGvtzOMeK6IOAHBINh5CrKrbc/RA0NNaa5c7kueTXDu12zVJLmxePACAPdaSNne47fBtNJJZVbck+eEkz2ytfXLqrbuT3FZVj6iq65PckOSt2xcTAIB9ssoURr+Y5KlJrqiq80lemKOnyR+R5DeqKkne3Fr72621d1XVq5O8O0dh9Oe21i72VXi2t8vl77aZhmPR3667LCDdmVf3N191enQ5gkP5nrvODWV17lccqqWdzNbas+Zsftkx+78oyYu2KRQAwKG4NPeRlcNnxR8AADpn7iAAgJ60jHfFH51M5uZt9ZXPJffocM07t/ICT950buzY8mH3kXsih0y4HACAzhnJBADoTVnxh/G6HFabN/XMPoRyFpXz5qtO//kPu6Hul+vrGpuu+3nXOPvFtcQ+MpIJANAjK/4AAEBHdDIBAOicTiYPMjvtyfSygEM3W859Kfchkgu7nr7yJeedh3nXiWvl5K1b787RfmutdvqzKzqZAAB0zoM/AAA9aW28K/4YyQQAoHNGMlnJvixPJ2+JfXNSec/Lrt1NyjD0+8FQzZuTGA6RTiYAQI+s+AMAAB0xkslci8I55y7cO9gQ2bywkzDUyRtq+9gHfV1fxx1zdvnJdexLGs1QzN5X4dDpZAIA9MiykgAA0BEjmQAAPRrrPJk6mSy1KM9xCDlFs2VYVKbZ73B5Pzmb3RlCe9hX0+1xl9fXptfFkHO1d809hjETLgcAoHNGMgEAetJSwuVjJnS6mb5DZKuEvqfLcNz0ILPTG83b1/ln1/Z59ZchpNEMIWy/6vnb1/MM69DJBADo0UhnMJKTCQBA93QyAQDonHB55ucSyZdZbpPcp1XyH9c97rxjzcvNWnR+nWuGaBf5hV1cC7N50n18h11es/K5WVsb7zyZRjIBAOickUwAgD6N9MkfI5kAAHTOSOaU6fyhfZ6v7qScdM7YJkvdXbbofK66DCXswknOPdllm+/z+ll27D7r7Li5eIGH0skEAOiRB38AAKAjRjJnCJOuZ9PQ1EmFmhZ9zrLPnw6vm9pqNUNYVvAQmBZnvk3SZfahLbq/jEPbgwd/quqWJC9NcirJz7fWXjxnn29N8mM5epTp3tbaf3fcMXUyAQBGrKpOJfmZJH89yfkkb6uqu1tr757a54YkL0jylNbax6vq8cuOK1wOADBuNyW5r7X2vtbap5PcleTWmX3+VpKfaa19PElaa/cvO6iRTACAnrQM4sGfK6rqnqnXZ1trZ6deX53kQ1Ovzyf5mpljfGmSVNX/naOQ+o+11n7tuA/VyaQTfU/tsWmu0qI8wW2mQ4KTdojLM65jiOXctkz7kC/KQfloa+3MMe/P6wXPZpI+LMkNSZ6a5Jok/2dVPbG19olFB9XJBADoS0uy+5HMZc4nuXbq9TVJLszZ582ttf+U5P1V9fs56nS+bdFB5WQCAIzb25LcUFXXV9XDk9yW5O6Zff5Vkq9Lkqq6Ikfh8/cdd1CdTACAEWutPZDkeUnOJXlPkle31t5VVT9eVc+c7HYuyceq6t1J3pDkB1trHzvuuMLldG6IS69tkz81xHywIRrKuT4Us0uhdjUP6T6156GV1XLDbGof5slsrb0uyetmtv3o1O8tyd+Z/KzESCYAAJ3TyQQAoHPC5fRqNrQ0b8m8VcJ/x+0jfMWh6nKJyV1dJ+uG+bss5/Rnd5XG08eyn+5hI7AH4fI+GMkEAKBzRjIBAHpTQ1jxZyeMZAIA0DkjmZyoeblHQ5zyaF2mNqFPm7Qv7fEz+qiLLq/5PvI8YQh0MgEA+uTBHwAA6IaRTAbl3IV79zZkDn1a5boYcrj1pK/tfUxhma6ffSs7x2jx4A8AAHRFJxMAgM4JlwMA9GmkD/7oZDI4i3KR+srn6mL6EPlT7MIm7W4fcxUP0bwlL6fvcc4Th0AnEwCgVx78AQCATuhkAgDQOeFy9sa8vKVNzTvG7Db5UOvp8vxwvHXzlqf332W7XtZGuizbqseat99JtOHpXHDXzAiM9MEfI5kAAHTOSCYAQJ9GOpK5tJNZVXcm+cYk97fWnjjZ9pgkr0pyXZIPJPnW1trHq6qSvDTJNyT5ZJLvbK39dj9Fh/WsG5Jad3/h9SPCf/1Tv7u3zfVueiLGYpVw+cuT3DKz7flJXt9auyHJ6yevk+QZSW6Y/NyR5Ge7KSYAAPtkaSeztfamJH88s/nWJK+Y/P6KJN80tf2V7cibkzy6qq7sqrAAAHulJWm1258d2fTBny9qrX04SSb/Pn6y/eokH5ra7/xk20NU1R1VdU9V3fORj13csBgAAAxR1w/+zOsuz013ba2dTXI2Sc6cfuRIU2LZxGzO37xlIbtYKpLNmc5oGPat/Q8lV3HevWR2+ybWvS+5jx2ONtJezqYjmX90OQw++ff+yfbzSa6d2u+aJBc2Lx4AAPto007m3Ulun/x+e5LXTm3/jjry5CR/cjmsDgDAeKwyhdEvJnlqkiuq6nySFyZ5cZJXV9VzknwwybdMdn9djqYvui9HUxg9u4cyAwDsj5GGy5d2Mltrz1rw1tPm7NuSPHfbQsEy83L+Vlkqsk/ypx7KnJnMM9supq/nk76OVvm8LspyXL6pHGYOlWUlAQDonGUlAQD6tMO5KndJJ5O9tmyqkV2En4YyDQvjNvQ2uMvyrXN/6Kqcxx1HmJxDpZMJANCjGumDP3IyAQDonE4mAACdEy7nYAwpB810Rp+xKG8Wdm16KqUul4+cts5xTfl1oFpGO0+mkUwAADpnJBMAoDdlCiM4RLteSeO4zx1jKL2r87Go7oQa97td9VX22Xax7HP6moZs0TRr05+1aDvsI+FyAAA6ZyQTAKBPHvwBAIBuGMnkoA05R2/M0xwtmqpl27owBQzzzLaLvq+9ZW1wXg7mcftwAIxkAgBAN3QyAQDonHA5AECfRhou18nkoO16nsxVjHU+vL6+87zjDvn8d2mM7WhV8+4F09deV23E3LjwGTqZAAB9aRntij9yMgEA6JyRTA7WWEKkLHfoUxsJwy636Pwv275K3a7atraZOmmsaTXsN51MAIAe1Ugf/BEuBwCgc0YyAQD6NNKRTJ1MGIAxLzF5Eg41H1N7Wd30VEVdTFvU5ZRHi6bd2ocp2OA4wuUAAHROJxMAgM7pZAIA0Dk5mRwkOUwcMrmYmxtK3c0ubzlrupyHPs8rh0snEwCgR+bJBACAjhjJZFCmQ0KbLr22z2bLP5TQHsOhTXTnuDB0H/U875izZVg0ndFxr4/7Wwai1a5LsBNGMgEA6JxOJgAAnRMuBwDoS4tlJWEI5BQ92LY5qhwO5797x+VwL7v21p1WaPoYx01fNO9zZz9/0WcftxSl9sMu6GQCAPRppCOZcjIBAOickUzYE6Yt2dy+r5jiHPdj1XaxbEWeLsuzSjnOXbh3YVh82d+u83nHHUObZBU6mQAAPbLiDwAAdMRIJgBAn0Y6kqmTCXtOntRq1s1hYxwWTS20bN91/3YTi445PVXRbF7pOuVYd5qj2e/rfsMywuUAAHTOSCYAQJ9GGi43kgkAQOeMZMIBkSe13D7Nmelcnqxt8nY3PVfT7XH2+l31mF3lhk6GGSXaAAAZh0lEQVTPw7lsn1X352j6IlMYAQBAR3QyAQDonHA5HBhTGsF29imlIll8zW/6HabD9scdY95nuecs0GrXJdgJI5kAAHROJxMAgM4JlwMA9GmkT5frZMIBkycFm5mXb7hqjuM6yzOusv+qn7domchNcjPn3TtWmSLNNGpM08kEAOiReTIBAKAjOpkAAHROuJyDsk0OEuOhnbCOddtLV7nQxx1n3nt950Ied3xzZi4hXA4AAN0wkgkA0Jc23gd/dDI5SPu2LBwwfIvCv4vuNV3dg+ZNC9R1KHrRcpTbTkk0HToXRh8f4XIAADpnJBMAoE8jDZcbyQQAoHNbjWRW1fcn+R9y1Ef/3STPTnJlkruSPCbJbyf59tbap7csJ7ABuU/Qv3Wvs+mcx9m/Pe69Ta37GbM57euWad60SqPPkTeSuZ6qujrJ9yQ501p7YpJTSW5L8pNJXtJauyHJx5M8p4uCAgCwP7YNlz8syWdX1cOSPCrJh5N8fZLXTN5/RZJv2vIzAADYMxuHy1trf1hV/yDJB5P8f0l+Pcnbk3yitfbAZLfzSa6e9/dVdUeSO5LkCVd7/ojujX0aI6Hy5YbeRkz5cphWmQqpz3O+Sgi76zD37BRJY2vTY50nc5tw+RcmuTXJ9UmuSvI5SZ4xZ9e5VdtaO9taO9NaO/O4x57atBgAAAzQNuHypyd5f2vtI621/5TkV5L81SSPnoTPk+SaJBe2LCMAAHtmm07mB5M8uaoeVVWV5GlJ3p3kDUm+ebLP7Uleu10RAQDYN9vkZL6lql6To2mKHkjyO0nOJvnXSe6qqp+YbHtZFwUFVjO2XCc4FKteu4vyJBctDXnccZfl/fZxP3GPGo+tnrhprb0wyQtnNr8vyU3bHBcA4GB48AcAALph7iAAgL608U5hpJPJQbOkGcc5xHZhbs3xOsT2zH4TLgcAoHNGMgEA+jTScLmRTEZhLKHDsXxP5psOlwqdssjNV53WPjgROpkAAHROuBwAoE/C5QAA0A0jmcDoyEdj6E66jd581Wk53T2pjHeeTCOZAAB0TicTAIDOCZczGlb/4RBpz4dll+dzCKtFTX//gwrfC5cDAEA3jGQCAPSlefAHAAA6YyST0Tl34V55bOyFRTlp2u9h6uu87lOeY1fl+8x3fm8nx2MzOpkAAH0SLgcAgG4YyQQA6NNIRzJ1Mhml6byfQ8lvG3quFevZ9nxqD/vjJO9Bq3zWPi8xOTvX56krd1ma/VJVtyR5aZJTSX6+tfbiBft9c5JfSvJXWmv3HHdM4XIAgBGrqlNJfibJM5LcmORZVXXjnP0+L8n3JHnLKsfVyQQA6FG13f6s4KYk97XW3tda+3SSu5LcOme//yXJTyX5s1UOKlzO6FlukqHY1xAlMHhXVNV0aPtsa+3s1Ourk3xo6vX5JF8zfYCq+uok17bWfrWqfmCVD9XJBADo0+4f/Ploa+3MMe/XnG1/Xuqq+gtJXpLkO9f5UOFyAIBxO5/k2qnX1yS5MPX685I8Mckbq+oDSZ6c5O6qOq7jqpMJADByb0tyQ1VdX1UPT3Jbkrsvv9la+5PW2hWttetaa9cleXOSZy57uly4HCYsNzkeQ8vDlYs5XkNpg7NmpwLaF4Msb8sQwuXHaq09UFXPS3IuR1MY3dlae1dV/XiSe1prdx9/hPl0MgEARq619rokr5vZ9qML9n3qKsfUyQQA6NGK0wgdHDmZAAB0zkgm7LlB5iCxsq7Pn/YwLEPNuVzVdPn7aFt9H5/d0skEAOiTcDkAAHRDJxMYrSGE59YNp+57+HVMnKvlpq9B9XV4hMsBAHrk6XIAAOiIkUwAgD6NdCRTJxP22BByCoEHO9TcwpuvOt3LPcd97HAJlwMA0DkjmQAAfWkRLoexG3qIS0iJ5KgdLGqr2gh9u9z2tDVWoZMJANCTmvyMkZxMAAA6p5MJAEDnhMthwOQ99e9yHQ89J5f9MJZ21Nd0RgdrpA/+GMkEAKBzRjIBAHpk7XIAAOiIkUwYGHlO49DledZm2AVzZrKMTiYAQJ+EywEAoBtGMgH2jPAk7BkjmQAA0A2dTAAAOidcDgDQl2aeTABOiJzKw+XcwmfoZAIA0DnhcgCAPgmXAwBAN3QyASKXDuhPtd3+7IpOJgAAndPJBACgc1s9+FNVj07y80memKO01r+Z5PeTvCrJdUk+kORbW2sf36qUMBJCtrt1uf5vvup0r8fnsPXdjthDHvzZyEuT/Fpr7cuSnE7yniTPT/L61toNSV4/eQ0AwIhsPJJZVZ+f5GuTfGeStNY+neTTVXVrkqdOdntFkjcm+eFtCgkAsK+s+LO+L0nykST/rKp+p6p+vqo+J8kXtdY+nCSTfx8/74+r6o6quqeq7vnIxy5uUQwAAIZmm5zMhyV5UpLvbq29papemjVC4621s0nOJsmZ048caR+fIRhC3pRcvcPnHANjs81I5vkk51trb5m8fk2OOp1/VFVXJsnk3/u3KyIAwJ5qA/jZkY07ma21/5DkQ1X1n082PS3Ju5PcneT2ybbbk7x2qxICALB3tl27/LuT/EJVPTzJ+5I8O0cd11dX1XOSfDDJt2z5GXDQhFGH59yFeztLo3B+OVTa9hpGmhS4VSeztfaOJGfmvPW0bY4LAMB+s+IPAACd2zZcDgDAApXxzpOpkwk7Ip/p8DnH43ZSub1DmIYN5tHJBADo00hHMuVkAgDQOZ1MAAA6J1wO0DG5mCTb5Uqu04am9z2p/ExtfD3VxhkvN5IJAEDnjGQCAPRlx+uH75JOJuyAUNPwXT5Hi8KPziHb0H42d/NVp9XfnhAuBwCgc0YyAQB6NNYVf4xkAgDQOSOZjFrf033IG9p/ziGb2kXb6Xs6o6FcD5e/21DKw3w6mQAAfRIuBwCAbhjJBADokQd/gM6d1BJvALPOXbj3IT+HYBfLaLIZnUwAADonXA4A0KeRhst1MhklIRaA9Qwp3D697Ov0dEamNhoWnUwAgL40D/4AAEBndDIBAOiccDn0SF4QMCTTeYur7r8vBp1rL1wOAADdMJIJANCTyngf/NHJZJSmp7/o8ngAQzfvfnXzVaf36j626N69T99hDITLAQDonJFMAIA+tXHGy41kAgDQOSOZjJJcTIDPWLQk45ByNafLNy+vfijlnGesD/4YyQQAoHM6mQAAdE64HACgLy2jXfFHJxPWNOS8H4BtDXp5xjy0fO7Jw6WTCQDQo7q06xLshpxMAAA6ZyQT1jRvmg8A+iNEvp90MgEA+jTSB3+EywEA6JxOJgAAnRMuhzXJBQIO0dCWahz6VErrsKwkAAB0xEgmAEBfWpI2zqFMI5kAAHTOSCajs02ej3xM4NAN4T437z49hHKxHp1MAIAeefAHAAA6YiQTViBMA9C/RelMe38PNpIJAADd0MkEAKBzwuUAAD2pjPfBH51MRmPTqYv2PhcIYI+5B+8vnUwAgL60ZsUfAADoipFMRuNyyGXdsPn0/sI2AP1xjz0sOpkAAD0a64M/wuUAAHTOSCYAQJ9GOpKpk8lobDqFUSJPCADWJVwOAEDnjGQCAPTIgz8AANCRrUcyq+pUknuS/GFr7Rur6vokdyV5TJLfTvLtrbVPb/s5sCvyMWE588kCs7oYyfzeJO+Zev2TSV7SWrshyceTPKeDzwAA2D8tyaW2258d2aqTWVXXJPmvk/z85HUl+fokr5ns8ook37TNZwAAsH+2DZf/r0l+KMnnTV4/NsknWmsPTF6fT3L1vD+sqjuS3JEkT7ja80f0b91lJYX8YHWuFziGB3/WU1XfmOT+1trbpzfP2XVu1bbWzrbWzrTWzjzusac2LQYAAAO0zRDiU5I8s6q+Ickjk3x+jkY2H11VD5uMZl6T5ML2xQQAYJ9sPJLZWntBa+2a1tp1SW5L8luttf8+yRuSfPNkt9uTvHbrUgIA7Klqu/3ZlT6SIX84yV1V9RNJfifJy3r4DBiURXme8tQAGKtOOpmttTcmeePk9/cluamL4wIA7L02zid/rPgDAEDndDIBAOicCSqhA+cu3GtZPQDm2uXDN7tkJBMAgM4ZyQQA6EuLFX9gLFYNZd981emVl6C8fNzLPwCwT6rqlqr6/aq6r6qeP+f9v1NV766qd1bV66vqi5cdUycTAGDEqupUkp9J8owkNyZ5VlXdOLPb7yQ501r7yiSvSfJTy44rXA4A0JNKUsOfJ/OmJPdN5jpPVd2V5NYk7768Q2vtDVP7vznJty07qJFMAIDDdkVV3TP1c8fM+1cn+dDU6/OTbYs8J8m/WfahRjIZpct5k6vkXO7b1ET7Vl6Ag3dp1wXIR1trZ455v+Zsmzv8WlXfluRMkr+27EN1MgEAxu18kmunXl+T5MLsTlX19CQ/kuSvtdY+teygwuUAAOP2tiQ3VNX1VfXwJLcluXt6h6r66iT/R5JnttbuX+WgRjIZtdlw8qLw+T6FnfeprABjMPQHf1prD1TV85KcS3IqyZ2ttXdV1Y8nuae1dneSn07yuUl+qaqS5IOttWced1ydTACAkWutvS7J62a2/ejU709f95g6mQAAfbHiDwAAdMdIJkyRzwgA3dDJBADoTUsG/uBPX4TLAQDonE4mAACdEy4HAOhRjTNabiQTAIDuGckEAOiTB38AAKAbOpkAAHROuBwAoC8tqUu7LsRuGMkEAKBzRjIBAPrkwR8AAOiGTiYAAJ0TLgcA6NM4o+VGMgEA6J6RTACAHpUHfwAAoBs6mQAAdE64HACgT8LlAADQDSOZAAB9aUlGuna5TiZ04OarTj/o9bkL9+6oJAAwDMLlAAB0zkgmAEBPKm2082TqZEIHzl2490Eh89nw+eV9AGAsdDIBAPo00pFMOZkAAHROJxMAgM4Jl8Oa5uVbLiIPEwDhcgAA6IhOJgAAnRMuBwDoi2UlgeOskocp/5J9cPNVp7XVKdPX9uV62aSO5h0Hxk4nEwCgR2Nd8UdOJgAAnTOSCQusM1URDM1x7XcfQ7uXyzy7hGsfnzH7+zbHWad+j/vM6eN09f335dyzv3QyAQD6JFwOAADdMJIJANCbNtqRTJ1MmEM+Jvugi3badVvfNs9vWXn27drsqrx9fO9Vc0BhU8LlAAB0zkgmAEBfWkYbLjeSCQBA54xkwsT0UnJ9zEkHl606j+IY2t4YvuM+Wve8yOFcYqRrlxvJBACgczqZAAB0buNweVVdm+SVSf6zHA0En22tvbSqHpPkVUmuS/KBJN/aWvv49kWF7s2GhDYN3QkV9W96WcGT/tw+P3O2zV3+rEMLIx/a9+HBVjm/Y75Plgd/1vZAkr/bWvvyJE9O8tyqujHJ85O8vrV2Q5LXT14DADAiG49kttY+nOTDk9//tKrek+TqJLcmeepkt1ckeWOSH96qlAAA+8pI5uaq6rokX53kLUm+aNIBvdwRffyCv7mjqu6pqns+8rGLXRQDAICB2HoKo6r63CS/nOT7Wmv/sapW+rvW2tkkZ5PkzOlHjrOLz4nqMyes77y9Q7doSp/pPMxF9XsSuX7LPuPchXsHvXzgrhzSd2F767QH99PDsFUns6o+K0cdzF9orf3KZPMfVdWVrbUPV9WVSe7ftpAAAHupJbk0zrG0jcPldTRk+bIk72mt/aOpt+5Ocvvk99uTvHbz4gEAsI+2Gcl8SpJvT/K7VfWOyba/l+TFSV5dVc9J8sEk37JdEdlH02HOeSESoZDxWTVUNm+/oYddh16+k6Qu6MKiqb32Uxvtgz/bPF3+fyVZlID5tE2PCwDA/rPiDwAAndv66XIAAI4hXA6bW5SHte72Vayam3PSuWG7WvZwl+TfAbCIcDkAAJ0zkgkA0KeRhsuNZAIA0DkjmXtujMt0DTUP8FDqd5mh1j9wmPb+3mrFHwAA6I5OJgAAnRMuH4CTCj92MW2QUOl4ONfASdr7sPhCLWmXdl2InTCSCQBA54xkAgD0yRRGAADQDSOZWxpL3tpYvuc2FtXRPuUZOc/ASdqn+yPr08kEAOiLeTIBAKA7RjIBAPo00gd/dDLXJGeNda3SZk46L0k7BqBvwuUAAHTOSCYAQJ+Ey3fnD975qJ19trAhQ7CsHW4aTte+gSG7+arTPS9b/N4ejsmqBtHJBAA4TG20I5lyMgEA6JxOJgAAnRtMuFzuGCzm+gAO1cHf31qSS5d2XYqdMJIJAEDndDIBAOjcYMLlAAAHydPlAADQDSOZAAB9MpIJAADd0MkEAKBzwuUAAL1pySXhcgAA6ISRTACAvrSkNSv+AABAJ3QyAQDonHA5AECfPPgDAADdMJIJANAnK/4AAEA3dDIBAOiccDkAQF9aSy6ZJxMAADphJBMAoE8e/AEAgG7oZAIA0DnhcgCAHjUP/gAAQDeMZAIA9KZ58AcAALqikwkAQOeEywEA+tKSXBIuBwCATuhkAgDQOeFyAIA+NfNkAgBAJ4xkAgD0pCVpHvwBAIBu6GQCANA54XIAgL605sEfAADoipFMAIAeefAHAAA60lsns6puqarfr6r7qur5fX0OAADbWdZvq6pHVNWrJu+/paquW3bMXjqZVXUqyc8keUaSG5M8q6pu7OOzAAAGrV3a7c8SK/bbnpPk4621v5TkJUl+ctlx+xrJvCnJfa2197XWPp3kriS39vRZAABsbpV+261JXjH5/TVJnlZVddxB+3rw5+okH5p6fT7J10zvUFV3JLlj8vJTv9le83s9leXQXJHko7suxB5QT6tRT6tTV6tRT6tTV6vZpp6+uMuCbOJP8/Fzv9lec8WOi/HIqrpn6vXZ1trZqddL+23T+7TWHqiqP0ny2BxzbvrqZM7r2T7o0arJlzubJFV1T2vtTE9lOSjqajXqaTXqaXXqajXqaXXqajX7Xk+ttVt2XYYVLO23rbjPg/QVLj+f5Nqp19ckudDTZwEAsLlV+m1/vk9VPSzJFyT54+MO2lcn821Jbqiq66vq4UluS3J3T58FAMDmVum33Z3k9snv35zkt1prx45k9hIun8Tqn5fkXJJTSe5srb3rmD85e8x7PJi6Wo16Wo16Wp26Wo16Wp26Wo166tmifltV/XiSe1prdyd5WZJ/XlX35WgE87Zlx60lnVAAAFibFX8AAOicTiYAAJ3beSfT8pOLVdUHqup3q+odl+e3qqrHVNVvVNV7J/9+4a7LuQtVdWdV3V9Vvze1bW7d1JF/PGlj76yqJ+2u5CdrQT39WFX94aRdvaOqvmHqvRdM6un3q+rm3ZT65FXVtVX1hqp6T1W9q6q+d7Jdm5pyTD1pUzOq6pFV9daqundSV39/sv36yZJ8750s0ffwyfa1l+w7BMfU08ur6v1TbeqrJttHee3trdbazn5ylFz675J8SZKHJ7k3yY27LNOQfpJ8IMkVM9t+KsnzJ78/P8lP7rqcO6qbr03ypCS/t6xuknxDkn+Tozm+npzkLbsu/47r6ceS/MCcfW+cXIOPSHL95No8tevvcEL1dGWSJ01+/7wkfzCpD21qtXrSph763SvJ505+/6wkb5m0lVcnuW2y/eeSfNfk9/8pyc9Nfr8tyat2/R12XE8vT/LNc/Yf5bW3rz+7Hsm0/OT6ppd1ekWSb9phWXamtfamPHR+rkV1c2uSV7Yjb07y6Kq68mRKulsL6mmRW5Pc1Vr7VGvt/Unuy9E1evBaax9urf325Pc/TfKeHK1uoU1NOaaeFhlzm2qttf938vKzJj8tydfnaEm+5KFtaq0l+w7BMfW0yCivvX21607mvGWMjrthjU1L8utV9fY6WoYzSb6otfbh5OiGn+TxOyvd8CyqG+3soZ43CTXdOZVyoZ6STMKUX52jERVtaoGZekq0qYeoqlNV9Y4k9yf5jRyN5H6itfbAZJfp+njQkn1JLi/Zd/Bm66m1drlNvWjSpl5SVY+YbBt1m9o3u+5krr1E0cg8pbX2pCTPSPLcqvraXRdoT2lnD/azSf5ikq9K8uEk/3CyffT1VFWfm+SXk3xfa+0/HrfrnG2jqas59aRNzdFau9ha+6ocrZ5yU5Ivn7fb5N/R1tVsPVXVE5O8IMmXJfkrSR6T5Icnu4+2nvbRrjuZlp88RmvtwuTf+5P8yxzdpP7ocmhg8u/9uyvh4CyqG+1sSmvtjyY39UtJ/mk+E74cdT1V1WflqOP0C621X5ls1qZmzKsnbep4rbVPJHljjnIIH11HS/IlD66PtZfsOzRT9XTLJDWjtdY+leSfRZvaS7vuZFp+coGq+pyq+rzLvyf5G0l+Lw9e1un2JK/dTQkHaVHd3J3kOyZPJT45yZ9cDoGO0Uz+0n+bo3aVHNXTbZOnXK9PckOSt550+XZhkvv2siTvaa39o6m3tKkpi+pJm3qoqnpcVT168vtnJ3l6jnJY35CjJfmSh7aptZbsOwQL6un/mfqfu8pR3up0mxrdtbevellWclVt/eUnx+SLkvzLSd73w5L8i9bar1XV25K8uqqek+SDSb5lh2Xcmar6xSRPTXJFVZ1P8sIkL878unldjp5IvC/JJ5M8+8QLvCML6umpk+lAWo5mMPgfk6QdLSH26iTvTvJAkue21i7uotw78JQk357kdye5YUny96JNzVpUT8/Sph7iyiSvqKpTORrQeXVr7Ver6t1J7qqqn0jyOznqtCcbLNl3IBbV029V1eNyFB5/R5K/Pdl/rNfeXrKsJAAAndt1uBwAgAOkkwkAQOd0MgEA6JxOJgAAndPJBACgczqZAAB0TicTAIDO/f89yfivf9SF9QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Quick check on the mask\n", "%pylab inline\n", "import matplotlib.pylab as plt\n", "plt.figure(figsize=[12,12])\n", "plt.pcolormesh(mask[0,:,:]) ; plt.colorbar()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "#---------------------------------------------------------------\n", "# sanity checks\n", "# compute water volume between 50C and 55C (never observed)\n", "wmass_testimpossible = lib_watermass.volume_watermass_from_ts(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),\\\n", " woa.s_an.transpose(),\\\n", " 50.,55.,0.,40.)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---Control checks---\n", "Impossible water mass, volume = 0.0\n", "All possible water, volume = 1.376204609569534e+18\n", "volume from scale factors = 1.3762046318318395e+18\n" ] } ], "source": [ "# compute all possible volume\n", "wmass_testall = lib_watermass.volume_watermass_from_ts(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),\\\n", " woa.s_an.transpose(),\\\n", " -10.,100.,0.,50.)\n", "\n", "# and compare to volume from metrics\n", "wmass_from_metrics = (dx * dy * mask * dz).sum()\n", "\n", "print('---Control checks---')\n", "print('Impossible water mass, volume =', wmass_testimpossible)\n", "print('All possible water, volume =', wmass_testall)\n", "print('volume from scale factors = ', wmass_from_metrics)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---Volume of ocean in temperature ranges---\n", "The percentage of ocean waters 20C < T < 40C is 1.6994850127711953 %\n", "The percentage of ocean waters 10C < T < 20C is 6.101128407070131 %\n", "The percentage of ocean waters 0C < T < 10C is 88.6711761110485 %\n" ] } ], "source": [ "#---------------------------------------------------------------\n", "# now let's learn some stuff\n", "# compute volume of water 20C < T < 40C\n", "wmass_20to40C = lib_watermass.volume_watermass_from_ts(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),\\\n", " woa.s_an.transpose(),\\\n", " 20.,40.,0.,40.)\n", "# compute volume of water 10C < T < 20C\n", "wmass_10to20C = lib_watermass.volume_watermass_from_ts(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),woa.s_an.transpose(),\\\n", " 10.,20.,0.,40.)\n", "# compute volume of water 0C < T < 10C\n", "wmass_0to10C = lib_watermass.volume_watermass_from_ts(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),woa.s_an.transpose(),\\\n", " 0.,10.,0.,40.)\n", "\n", "print('---Volume of ocean in temperature ranges---')\n", "print('The percentage of ocean waters 20C < T < 40C is ', 100 * wmass_20to40C / wmass_from_metrics, '%')\n", "print('The percentage of ocean waters 10C < T < 20C is ', 100 * wmass_10to20C / wmass_from_metrics, '%')\n", "print('The percentage of ocean waters 0C < T < 10C is ', 100 * wmass_0to10C / wmass_from_metrics, '%')" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Solution from v1 is 1.6994850127711953 %\n", "Solution from v2 is 1.6994850127711953 %\n", "Solution from v2 without nx,ny,nz is 1.6994850127711953 %\n", "Solution from v3 without nx,ny,nz is 1.6994850127711953 %\n" ] } ], "source": [ "# note on I/O and dimensions\n", "\n", "#---------------------------------------------------------------\n", "\n", "wmass_20to40C = lib_watermass.volume_watermass_from_ts(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),\\\n", " woa.s_an.transpose(),\\\n", " 20.,40.,0.,40.)\n", "\n", "print('Solution from v1 is ', 100 * wmass_20to40C / wmass_from_metrics, '%')\n", "\n", "wmass_20to40C = lib_watermass.volume_watermass_from_ts_v2(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),\\\n", " woa.s_an.transpose(),\\\n", " 20.,40.,0.,40.,nx,ny,nz)\n", "\n", "\n", "print('Solution from v2 is ', 100 * wmass_20to40C / wmass_from_metrics, '%')\n", "\n", "wmass_20to40C = lib_watermass.volume_watermass_from_ts_v2(dx.transpose(),dy.transpose(),dz.transpose(),\\\n", " woa.t_an.transpose(),\\\n", " woa.s_an.transpose(),\\\n", " 20.,40.,0.,40.)\n", "\n", "print('Solution from v2 without nx,ny,nz is ', 100 * wmass_20to40C / wmass_from_metrics, '%')\n", "\n", "wmass_20to40C = lib_watermass.volume_watermass_from_ts_v3(dx,dy,dz,\\\n", " woa.t_an,\\\n", " woa.s_an,\\\n", " 20.,40.,0.,40.)\n", "\n", "print('Solution from v3 without nx,ny,nz is ', 100 * wmass_20to40C / wmass_from_metrics, '%')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "gist_info": { "gist_id": null, "gist_url": null }, "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.5" } }, "nbformat": 4, "nbformat_minor": 2 }