0-unit_conversions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0b8477ac-b117-458a-8542-e209d5f499b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "table {align:left;display:block}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<style>\n",
    "table {align:left;display:block}\n",
    "</style>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "323de01b-930e-443b-8c6e-dd8adf7399e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79c8ded0-f0c2-4961-aa3f-7079adf4481e",
   "metadata": {},
   "source": [
    "# Units\n",
    "\n",
    "Units of measurements are extremely important, especially while comparing with the results of others.\n",
    "\n",
    "This notebook is devoted to unit conversions.\n",
    "\n",
    "\n",
    "# Productivity \n",
    "\n",
    "<img align=\"right\" src=\"productivity_parts.png\" width=300 height=300/>\n",
    "\n",
    "\n",
    "\n",
    "| quantile | meaning | units |\n",
    "| --- | --- | --- |\n",
    "| GPP | Gross primary productivity | $\\mu$mol CO$_2$ m$^{-2}$ s$^{-1}$ |\n",
    "| AR | Autotrophic (plant) respiration | $\\mu$mol CO$_2$ m$^{-2}$ s$^{-1}$ |\n",
    "| NPP | Net primary productivity | $\\mu$mol CO$_2$ m$^{-2}$ s$^{-1}$ |\n",
    "| DMP | Dry matter productivity | kg DM ha$^{-1}$ day$^{-1}$ |\n",
    "| AGBP | Aboveground biomass productivity | kg DM ha$^{-1}$ day$^{-1}$ |\n",
    "| yield | crop yield | kg FW ha$^{-1}$ season$^{-1}$ |\n",
    "\n",
    "\n",
    "DM - dry matter – CH$_2$O or C$_6$H$_{12}$O$_6$\n",
    "\n",
    "FW – fresh matter (weight)\n",
    "\n",
    "Image from Blatchford, Megan L, Chris M Mannaerts, Yijian Zeng, Hamideh Nouri, and Poolad Karimi. 2019. “Status of Accuracy in Remotely Sensed and In-Situ Agricultural Water Productivity Estimates: A Review.” Remote Sensing of Environment 234 (September 2018): 111413. https://doi.org/10.1016/j.rse.2019.111413\n",
    "\n",
    "## GPP to NPP\n",
    "\n",
    "NPP = GPP - AR\n",
    "\n",
    "AR is usually not known and a constant fraction is used (NPP = 0.5 * GPP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7efdc16-d934-4a1d-94b2-b7e33dff4ef1",
   "metadata": {},
   "source": [
    "## GPP, NPP, AR\n",
    "\n",
    "### umol to g C\n",
    "\n",
    "here chemistry has to be remembered \n",
    "\n",
    "1. mol is the quantity of matter that containts 6.02E23 particles\n",
    "2. molar mass (M) shows how many grams are in one mol ([periodic table](https://ptable.com/#Properties))\n",
    "3. the numbers in formulae show how many mols of certain element are there\n",
    "    - CO2 has 1 mol of C (M=12 g mol-1) and 2 mols of O (M = 16 g mol-1)\n",
    "    - M(CO2) = M(C) + M(O) * 2 = 12 + 16 * 2 = 44 g mol-1\n",
    "4. **note, we work with micromoles, 1e-6 mol**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "be692ae5-4dc8-4df8-9b23-4033166211a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpp_umol = 10 # umol CO2 m-2 s-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "34259bea-30e3-4df0-920c-7335677c0a6b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00011999999999999999"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# finish the equation\n",
    "gpp_gC = gpp_umol * 12 * 1e-6 # g C m-2 s-1\n",
    "gpp_gC"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f250ae7-4b6f-41c8-beab-b1081447d0df",
   "metadata": {},
   "source": [
    "### s-1 to day-1\n",
    "\n",
    "1 day = 24 hours = 24 * 60 minutes = 24 * 60 * 60 seconds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "4a85a784-a60f-4dcc-b35c-892766f45b06",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpp_gC_s = 0.00012  # g C m-2 s-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "5abab09e-7714-4be1-91c7-b7c36d01d00f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.368"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gpp_gC_day = gpp_gC_s * 24 * 60 * 60  # g C m-2 day-1\n",
    "gpp_gC_day"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75cb5ae2-5ff7-4143-9f73-4dad8d49d840",
   "metadata": {},
   "source": [
    "### m-2 to ha-1\n",
    "\n",
    "1 ha = 10000 m2 = 1e4 m2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "45119fe4-425b-4d85-b167-dc2ca522e592",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpp_gC_day_m2 = 10.368  # g C m-2 day-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "3d24c99d-2ea9-4388-86fd-929c565c522c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "103680.0"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gpp_gC_day_ha = gpp_gC_day_m2 * 1e4  # g C ha-1 day-1\n",
    "gpp_gC_day_ha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "79b6a2f2-7c0e-4c31-9dfb-3809268034e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "230.4"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gpp_gC_day_m2 / 0.045  # eq 4 in Blatchford et al., 2019 NPP2DMP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "999d09bc-1425-4a10-9358-13d731e1d310",
   "metadata": {},
   "source": [
    "## DMP\n",
    "\n",
    "### g C kg DM\n",
    "\n",
    "DM - dry matter – CH$_2$O (general carbohydrate) or C$_6$H$_{12}$O$_6$ (glucose)\n",
    "\n",
    "the mass to mass conversion is done through quantity (mols)\n",
    "\n",
    "1 kg = 1e3 g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d54fe55c-d750-4a68-8e66-b54e756c68f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpp_gC_day_ha = 103680  # g C ha-1 day-1\n",
    "gpp2npp = 0.5\n",
    "npp_gC_day_ha = gpp_gC_day_ha * gpp2npp  # g C ha-1 day-1\n",
    "npp_molC_day_ha = npp_gC_day_ha / 12 # mol C ha-1 day-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "193daf6e-b5cc-48c5-a808-253318fb4d9d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "129600.0"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M_CH2O = 12 + 2 + 16\n",
    "dm_molCH2O_day_ha = npp_molC_day_ha / 1 # mol CH2O ha-1 day-1\n",
    "dm_gCH2O_day_ha = dm_molCH2O_day_ha * M_CH2O\n",
    "dm_gCH2O_day_ha"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9919822-0fea-4f33-8789-225fb38ab9a3",
   "metadata": {},
   "source": [
    "what if we want to compute it for C$_6$H$_{12}$O$_6$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "aa78576f-c31f-488a-b945-d2901385db2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "129600.0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M_C6H12O6 = ... # g mol-1\n",
    "M_C6H12O6 = 12*6 + 12 + 16*6 # g mol-1\n",
    "\n",
    "dm_molC6H12O6_day_ha = npp_molC_day_ha / 6 # mol C6H12O6ha-1 day-1\n",
    "dm_gC6H12O6_day_ha = dm_molC6H12O6_day_ha * M_C6H12O6\n",
    "dm_gC6H12O6_day_ha"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3da36f16-8541-418c-be34-5133c8d7c996",
   "metadata": {},
   "source": [
    "## AGBP\n",
    "\n",
    "There are many studies that estimate above ground biomass from vegetation indices, RGB images and lidar, but we are talking about AGBP/DMP ratio. \n",
    "\n",
    "In other words, crop-specific shoot to total mass ratio, above ground biomass fraction (AGB fraction, f_agb)\n",
    "\n",
    "As you may expect, this is a dynamic value over the growing season.\n",
    "\n",
    "The example computations provided here are from \n",
    "\n",
    "Baret, F., A. Olioso, and J. L. Luciani. 1992. “Root Biomass Fraction as a Function of Growth Degree Days in Wheat.” Plant and Soil 140 (1): 137–44. https://doi.org/10.1007/BF00012815."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "8db8297b-6aef-4185-8a05-dee97c75b32e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'aboveground biomass partitioning fracton')"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxlklEQVR4nO3dd3wVZfb48c9JKAmBECCAEAgdBJQaAcHeXQt2QVaxreKua9v+c3dtu6uuZVe/6iqLKFbsil1xxUbvHaQTCCGUJKS38/tjJnqJ4WZCcu/c5J7363VfmZk7MznMDXPuzDzPeURVMcYYE71i/A7AGGOMvywRGGNMlLNEYIwxUc4SgTHGRDlLBMYYE+Wa+B1AbSUnJ2v37t39DsMYYxqURYsW7VHV9tW91+ASQffu3Vm4cKHfYRhjTIMiIlsP9Z7dGjLGmChnicAYY6KcJQJjjIlylgiMMSbKeXpYLCIpQLfA9VX161AFZYwxJnxqTAQi8iBwObAaKHcXK2CJwBhjGgEvVwQXAP1UtTjEsRhjjPGBl0SwCWgKWCIwxpgwKSotJyOniIzsQna6P4ekJnF8n2r7hNWJl0RQACwVkS8ISAaqeku9R2OMMVGgpKyCzNwidmYXsiu3iJ3ZRWTkFP7wMyOniH35JT/ZbtKJvXxLBDPclzHGGA+KSsvZkV3Itr0FbNvnvHbsL3RO9jlF7MkrpuqYYK3imtC5dTydkuIY1CWJzq3j6JQU/8PPTq3jiGsaG5J4a0wEqjpNRJoBfd1F61S1NCTRGGNMA6Cq7M0vYdu+ArbvK2Db3gK2uif87fsK2JVbdNCJPq5pDClJ8XRqHc9JHVsddIKv/NmyuX8Vf7y0GjoJmAZsAQToKiITrfmoMaYxq6hQMnKL2Lg7jy178w/6dr99XwH5JeUHrd8xsTmpbVtwbK92pLZtcdCrfavmiIhP/5KaeUlBjwBnqOo6ABHpC7wKDA9lYMYYEw5FpeVs3VvAht15bMxyXht257EpK5/C0h9P9nFNY344sVc92Xdp04L4ZqG5bRMOXhJB08okAKCq60WkqZedi8hZwGNALDBFVR+o8n5r4CUg1Y3lYVV9zmvwxhjjVXZBScDJPv+H6e37CqgIuI2TkhRP7w4tGdmjHb06JNCrfUt6JCfQIcK/1deFl0SwUESeBV505ycAi2raSERigSeB04F0YIGIzFDV1QGr/QpYrarniUh7YJ2IvKyqP31cbowxHhSUlLFu1wFWZ+SyJiOX9bucE/7egFY4zZrE0DM5gaNSWjN2SAq9O7SkV/sEeia3bNDf7A+Xl0RwE84J+xacZwRf45zgazIC2KCqmwBEZDowFqeHciUFWomTZlsC+4Ayz9EbY6KWqpKZW8yajFxWu681Gbls3pP/w4PaVnFNOPKIVpwxsCO92rf84ZXSJp7YmMb57f5weEkEk1T1UeDRygUicivOLZ9gUoDtAfPpwMgq6zyB0zR1J9AKuFxVK6ruSERuAG4ASE1N9RCyMaYxKS2vYGNWHqt35v5w4l+TceCgtvZd28YzoFMi5w/uzIBOiQzonEhKUnyjvZ1Tn7wkgon89KR/dTXLqqru6FdpOcuZwFLgFKAX8LmIfKOquQdtpDoZmAyQlpZWdR/GmEaktLyCtRkHWLJ9P8vTc1iTkcv3mXmUlDvfEZs1iaFfx1ac3r8jAzon0r9TIkd2akVinKdHl6Yah0wEIjIeuALoISKBHcpaAXs97Dsd6Bow3wXnm3+ga4AHVFWBDSKyGTgSmO9h/8aYRiAzt4gl2/azZFs2i7c5J//iMuek3y6hGQM6J3LNmO4/nPR7JifQJNYq6NenYFcEs4EMIBmnCWmlA8ByD/teAPQRkR7ADmAcTmIJtA04FfhGRDoC/XBqGxljGqGi0nJW7cxhybZs97WfnTlFADSLjWFgSiITRnZjaGoSQ1OT7NZOmBwyEajqVmCriEwAdqpqEYCIxON8u98SbMeqWiYiNwOf4jQfnaqqq0Rkkvv+08B9wPMisgLnVtIfVHVP3f9Zxhi/qSrp+wtZ7H7bX7JtP6szciktd+7udmkTz/Dubbm+q3PSH9A5keZNoq/FTiQQrVrwouoKIguB0ZVNOt1yE9+p6jFhiO8n0tLSdOHChX78amNMEKrKpj35zN64lzkb9zB/8z725DkPc+ObxjK4a2uGprZhaNckhqQm0aFVnM8RRxcRWaSqadW95+VhcZPAdv2qWuImA2NMlEvfX+Ce+Pcye+MeMnOdAsUpSfGc0Lc9w7u1YWjXNvTt2NLu60cwL4kgS0TOV9UZACIyFrDbN8ZEoawDxczZtJfZG/Ywe+Netu0rACC5ZTOO7ZXM6F7tGO2WX7B7+w2Hp34EwMsi8gTOffztwFUhjcoYExFyCkqZu/nHb/zrM/MAp6PWqJ7tuGZMd0b3SqZvx5Z24m/AvJSh3giMEpGWOM8UDoQ+LGOMH8rKK1iwZT+z1u9m9oa9rNyZg6pTcO2Y7m25cGgXxvRux8DOra1nbiPiqQC2iJwDDATiKrO+qt4bwriMMWFyoKiUr9ZnMXN1Jl+uyyKnsJSmscLQ1DbcemofRvdKZnDX1taipxHzMh7B00AL4GRgCnAJ1uHLmAYtfX8BX6zZzcw1mczdtJfScqVNi6ac1r8jpw/owPF92pPg40ApJry8fNKjVXWQiCxX1XtE5BHg7VAHZoypP6rKih05zFydyedrdrMmw6ni0rN9AteO6cFpAzoyLLWN3e6JUl4SQZH7s0BEOuOUl+gRupCMMfWhqLScOZv2MnN1JjPXZJKZW0yMQFq3tvy/nx3Jqf2dipzGeEkE74tIEvAQsBincNx/QxmUMebw5BSW8vnqTGauzuTr77MoKCmnRbNYTuzbntP6d+TkIzvQNsG6AZmDBU0EIhIDfKGq2cBbIvIBEKeqOeEIzhhTs+Kycmaty+LdJTv4Yu1uSsoqOCIxjouGpXBa/46M6tmOuKb2oNccWtBEoKoV7jOBY935YqA4HIEZYw5NVVm0dT9vL9nBh8szyCkspV1CM64YkcoFQ1MY3KW1tes3nnm5NfSZiFwMvK01FSYyxoTUht15vLd0B+8s2UH6/kLimsZw5sAjuGBoCsf1TqaplXEwh8FLIrgDSADKRKQIp3exqmpiSCMzxgBOWYf3l+3k3aU7WJ6eQ4zAmN7J3HF6X84YeAQtrZmnqaNgA9OMUdXvgPaVJaiNMeFRUFLGZ6syeWfJDr7dsIfyCuWolET+fE5/zh/cmQ6JVrnT1J9gXyUeB4bjDFAzLDzhGBO9yiuUbzfs4d0lO/h01S4KSspJSYpn0ok9uWBICn06tvI7RNNIBUsEpSLyHNBFRB6v+qaq3hK6sIyJHjkFpby+cDsvzt3Ktn0FJMY1YeyQFC4cmkJatzbEWCcvE2LBEsG5wGk4A8svCk84xkSPtbtymTZ7C+8s2UFRaQXHdG/D787sxxkDO1pdHxNWwYaq3ANMF5E1qrosjDEZ02iVlVfw2epMps3ewrzN+2jeJIYLhqRw1ehuDOzc2u/wTJTyUobakoAxdbQnr5jp87fx8rxtZOQU0aVNPH86+0guP6YrSS2sp6/xl7U7MyaElm3PZtrsLXywPIOS8gqO653MvWOP4pQjO1iBNxMxLBEYU8+Ky8r5aEUG02ZvZen2bBKaxTJuRFeuOrYbvTtYyx8TebyMR3BHNYtzgEWqurTeIzKmgdqVU8TL87by6vxt7MkroWdyAnefN4CLh3ehVVxTv8Mz5pC8XBGkua/33flzgAXAJBF5Q1X/GargjGkItu8r4KlZG3hjYTrlqpzSrwMTR3fnuN7J1vTTNAheEkE7YJiq5gGIyF3Am8AJOM1KLRGYqBSYAGJEGD8ilV8c35PUdi38Ds2YWvGSCFKBkoD5UqCbqhaKiFUiNVGnagK4YmQqN53Ui06t4/0OzZjD4iURvALMFZH33PnzgFdFJAFYHbLIjIkw6fsLePLLjbyxcLslANOoeOlHcJ+IfAyMwak8OklVF7pvTwhlcMZEAksAprHz2nx0CbCzcn0RSVXVbSGLypgIUJkA3ly0HcESgGm8vDQf/TVwF5AJlOOORwAMCm1oxvijagIYP8ISgGncvFwR3Ar0U9W9oQ7GGD+l7y/gqVnOLSBBGHeMkwA6J1kCMI2bl0SwHacDmTGN0p68Yv71+XpetwRgopSXRLAJmCUiHxIwcL2qPhqyqIwJg7LyCl6Ys5V/zVxPYUk540Z05Zcn9bYEYKKOl0SwzX01c1/GNHhzNu7l7hmrWJd5gOP7JHPXeQPp3aGl32EZ4wsvzUfvCUcgxoTDzuxC/v7RGj5cnkGXNvE8c+VwzhjQERErBWGiV7DB6/+tqreJyPs4rYQOoqrnhzQyY+pRUWk5U77ZxJNfbqRCldtP68uNJ/YkrqmNBGZMsCuCF92fD4cjEGNCQVX5Ys1u7v1gNdv2FXD2UUdw5zn96dLG6gEZUynYUJWL3J9fhS8cY+rP5j353PP+Kmaty6J3h5a8dN1IjuuT7HdYxkQcLx3KxgB3A93c9QVQVe3pYduzgMeAWGCKqj5QzTonAf8GmgJ7VPVEz9EbU4384jKe+HIDz36zmWZNYvjzOf2ZOLo7TWNj/A7NmIjkpdXQs8DtOCWny73uWERigSeB04F0YIGIzFDV1QHrJAFPAWep6jYR6VCL2I05iKoyY9lO7v9oLbtyi7h4WBf+cHY/OrSK8zs0YyKal0SQo6ofH8a+RwAbVHUTgIhMB8ZycMXSK4C3K+sWqeruw/g9xrB6Zy53v7+K+Zv3cXRKa56cMIzh3dr4HZYxDYKXRPCliDwEvM3BHcoW17BdCk6v5ErpwMgq6/QFmorILKAV8JiqvlB1RyJyA3ADQGpqqoeQTbQoKCnjn5+s44U5W2gd35T7Lzqay9K62sDwxtSCl0RQefJOC1imwCk1bFfd/8SqzVCbAMOBU4F4YI6IzFXV9QdtpDoZmAyQlpb2k6asJjot2bafO15fxpa9+fx8ZDd+c0ZfklpYn0djastLh7KTD3Pf6UDXgPkuOKWsq66zR1XzgXwR+RoYDKzHmEMoLa/g//63gSe/3MARiXG8cv0oju3Vzu+wjGmwvLQaao1ThvoEd9FXwL2qWlMhugVAHxHpAewAxuE8Ewj0HvCEiDTBKV8xEviX9/BNtNmYlccdry1lWXoOFw1L4e7zB5IY19TvsIxp0LzcGpoKrAQuc+evBJ4DLgq2kaqWicjNwKc4zUenquoqEZnkvv+0qq4RkU+A5UAFThPTlYf3TzGNmary0tyt/P2jNcQ1jeWpCcP42dGd/A7LmEZBVIPfcheRpao6pKZl4ZKWlqYLFy6seUXTaOzOLeJ3by7nq/VZnNC3PQ9dMoiOidYk1JjaEJFFqppW3XterggKReQ4Vf3W3dkYoLA+AzTmUD5ekcGf3llBUWk5940dyM9HdbMCccbUMy+J4CZgmvusQIB9wNWhDMqY3KJS7p6xircX72Bwl9Y8evkQerW3MtHGhIKXVkNLgcEikujO54Y6KBPd5m7ay29eX8au3CJuPbUPN5/S28pDGBNCwcpQ/1xVXxKRO6osB2yEMlP/isvKefSz9Uz+ZhPd2rbgzUnHMjTVegcbE2rBrggS3J+tqnnPOnWZerUmI5fbX1vK2l0HmDAylTvP6U+LZl7uXBpj6ipYGepn3MmZqvpd4HvuA2Nj6qyiQpny7SYe/nQ9ifFNmXp1Gqcc2dHvsIyJKl6+cv0fMMzDMmNqZV9+CTe/spjZG/dy5sCO/OPCo2nXsrnfYRkTdYI9IzgWGA20r/KcIBGng5gxh23trlyun7aQ3QeK+efFg7g0rYs1CzXGJ8GuCJoBLd11Ap8T5AKXhDIo07h9tmoXt7+2lITmTXj9xmMZ0jXJ75CMiWrBnhF8BXwlIs+r6tYwxmQaKVXlyS838PBn6xncpTWTr0qzHsLGRIBgt4b+raq34RSF+0krIVU9P5SBmcalsKSc3725jA+WZ3DBkM48cPEg4praHUZjIkGwW0Mvuj8fDkcgpvHamV3IDS8uZNXOXP549pHceEJPex5gTAQJdmtokTs5RFUfC3xPRG7FKUdtTFCLtu7nxhcXUVRazrMTrWmoMZHIS7/9idUsu7qe4zCN0BsLtzN+8lwSmsfyzi9HWxIwJkIFe0YwHmcgmZ4iMiPgrVbA3lAHZhqusvIK7v94Lc9+u5kxvdvx5BXDbAhJYyJYsGcEs4EMIBl4JGD5AZyBZIz5iZzCUn796hK+Xp/F1aO7c+c5/a1gnDERLtgzgq0ikg7ku01JjQlqY1Yev5i2kO37C7j/oqMZPyLV75CMMR4ELTGhquUiUiAirT2MUWyi2Ffrs7j5lcU0i43h5etHMaJHW79DMsZ45KXWUBGwQkQ+B/IrF6rqLSGLyjQYqsqz327mHx+toW/HVkyZmEaXNi38DssYUwteEsGH7suYgxSVlnPnOyt5a3E6Zw08gkcuG0xCcysdbUxD42WEsmnhCMQ0LAeKSrnu+YXM37KPW0/tw62n9iEmxjqJGdMQ1ZgIRKQPcD8wAPihMIyq9gxhXCaC5RSUctVz81m1I4fHxw/l/MGd/Q7JGFMHXtr1PQf8BygDTgZe4MfyEybK7Msv4Yopc1mzM5f//Hy4JQFjGgEviSBeVb8ARFW3qurdwCmhDctEoqwDxYyfPJcNu/OYfNVwTh9gPYWNaQw8tRoSkRjgexG5GdgBdAhtWCbS7Mop4oopc8nILuK5q49hdO9kv0MyxtQTL1cEtwEtgFuA4cCVVF9/yDRS6fsLuOyZOezOLeaF60ZYEjCmkfHSamgBgHtVcIuqHgh5VCZibN2bzxX/nUduUSkvXjeCoalt/A7JGFPParwiEJE0EVmBU19ohYgsE5HhoQ/N+G1jVh6XPzOX/JIyXv3FKEsCxjRSXp4RTAV+qarfAIjIcTgtiQaFMjDjr3W7DjBhyjxAmX7DKI48ItHvkIwxIeLlGcGByiQAoKrf4lQgNY3Uyh05jJs8hxiB6Tcca0nAmEbOyxXBfBF5BngVUOByYJaIDANQ1cUhjM+E2bLt2Vz57DxaNm/CK78YRffkBL9DMsaEmJdEMMT9eVeV5aNxEoP1KWgkFm7Zx9XPLaBNQlNeuX4UXdta8ThjooGXVkMnhyMQ4685G/dy3bQFHJEYx8u/GEmn1vF+h2SMCRMbOso4o4k9N5+UpHim3zjKkoAxUcZqBke5L9ZkctNLi+nVoSUvXTeCdi2b+x2SMSbM7Iogin28IoMbX1xE/06tePUXIy0JGBOlvHQou1REWrnTfxaRtytbDHnY9iwRWSciG0Tkj0HWO0ZEykXkEu+hm7r4YPlObn51CYO7JvHi9SNJatHM75CMMT7xckXwF1U94HYkOxOYhlOWOigRiQWeBM7GGctgvIgMOMR6DwKf1iZwc/hmb9zD7a8tZXhqG164dgSJcU39DskY4yMviaDc/XkO8B9VfQ/w8vVxBLBBVTepagkwHRhbzXq/Bt4CdnvYp6mj9ZkHuPHFRfRITuC/E9NsaEljjKdEsMPtUHYZ8JGINPe4XQqwPWA+3V32AxFJAS4EnvYWrqmLzNwirp46n/imsTx3zQhax9uVgDHG2wn9MpzbNmepajbQFvidh+2qG8BWq8z/G/iDqpZXs+6POxK5QUQWisjCrKwsD7/aVJVXXMa1zy8gp7CUqVcfQ0qSNRE1xji83BfoBHyoqsUichJOsbkXPGyXDnQNmO8C7KyyThowXUQAkoGfiUiZqr4buJKqTgYmA6SlpVVNJqYGZeUV/OrlxazddYApE9M4KqW13yEZYyKIlyuCt4ByEekNPAv0AF7xsN0CoI+I9BCRZsA4YEbgCqraQ1W7q2p34E2cKqfv1iJ+UwNV5c/vruSr9Vn87YKjOLmfDS5njDmYl0RQoaplwEXAv1X1dpyrhKDcbW7Gua20BnhdVVeJyCQRmVSXoI13T83ayPQF27n55N6MH5HqdzjGmAjk5dZQqYiMB64CznOXeXrKqKofAR9VWVbtg2FVvdrLPo137yxJ56FP13Hh0BR+c0Zfv8MxxkQoL1cE1wDHAn9X1c0i0gN4KbRhmbqavXEPv39zOcf2bMeDFw/CfQ5jjDE/4aX66Gqcgesr5zcDD4QyKFM3lX0FurdL4Okrh9OsiVUSMcYcWo2JQET6APfj9A6Oq1yuqj1DGJc5TIF9BZ6/1voKGGNq5uWr4nM4JSXKgJNxmo6+GMqgzOGp7CuQbX0FjDG14CURxKvqF4Co6lZVvRsblSziBPYVeHLCMOsrYIzxzEuroSIRiQG+F5GbgR2ANUaPIIF9Be6/6GjrK2CMqRUvVwS3AS1wHhgPB64EJoYwJlNL1lfAGFMXXloNLXAn83CakpoIYn0FjDF15aXVUBpwJ9AtcH1VHRTCuIwH1lfAGFMfvDwjeBmn2ugKoCK04RivrK+AMaa+eEkEWao6o+bVTLhYXwFjTH3ykgjuEpEpwBdAceVCVX07ZFGZQyotr+CmlxaRXVjK6zcea30FjDF15iURXAMciVNorvLWkAKWCHzwyGfrWbwtmyeuGGp9BYwx9cJLIhisqkeHPBJTo1nrdvP0VxsZPyKVcwd19jscY0wj4eUJ41wRGRDySExQmblF3PH6Mvp1bMVd59nHYYypP16uCI4DJorIZpxnBAKoNR8Nn/IK5ZZXl1BYUs6TE4YS1zTW75CMMY2Il0RwVsijMEE9/sX3zNu8j4cvHUzvDq38DscY08jUeGtIVbcCSTijk50HJLnLTBjM3riHx//3PRcNS+GS4V38DscY0wjVmAhE5FacTmUd3NdLIvLrUAdmYE9eMbdNX0qP5ATuG3uU3+EYYxopL7eGrgNGqmo+gIg8CMwB/i+UgUW7igrljteXkV1YyrRrR5DQ3MtHZYwxteel1ZAA5QHz5e4yE0LPfL2Jr9dncdd5A+jfKdHvcIwxjZiXr5nPAfNE5B13/gLg2ZBFZFi0dR8Pf7aOcwZ14gorK22MCTEvZagfFZFZOM1IBbhGVZeEOrBolV1Qwq9fWUJKUjz3X3S0VRQ1xoTcIROBiCSqaq6ItAW2uK/K99qq6r7QhxddVJXfvrGcrLxi3pw0msQ4KyZnjAm9YFcErwDnAotwagsFfjVVoGcI44pKz323hZlrMvnLuQMY3DXJ73CMMVHikIlAVc91f/YIXzjRa3l6Nvd/vIbT+nfk2jHd/Q7HGBNFPLVJFJGLcJ4RKPCNqr4byqCiTW5RKTe/soT2LZvz8KU20pgxJry8dCh7CpiEM0LZSmCSiDwZ6sCiharyp7dXsCO7kMfHDyWpRTO/QzLGRBkvVwQnAkepqgKIyDScpGDqwavzt/Ph8gx+f1Y/0rq39TscY0wU8tKhbB0Q2Ji9K7A8NOFElzUZudzz/iqO75PMpBN6+R2OMSZKBWs++j7OM4HWwBoRme/OjwRmhye8xiu/uIybX1lMYnxT/nX5EGJi7LmAMcYfwW4NPRy2KKLQX99bxaY9+bx83UiSWzb3OxxjTBQL1nz0q3AGEk3eWpTOW4vTueXUPozunex3OMaYKOflGYGpRxt25/GX91Yyskdbbj21j9/hGGOMJYJwqqhQfv/mMpo3ieGxcUOJtecCxpgIYIkgjF6ev43F27L58zkDOKJ1nN/hGGMMELzV0AqcVkLVssHra2d3bhH//HgtY3q346JhKX6HY4wxPwjWauhc9+ev3J8vuj8nAAVedi4iZwGPAbHAFFV9oMr7E4A/uLN5wE2quszLvhuae95fTXF5BX+7wEpLG2MiS7BWQ1sBRGSMqo4JeOuPIvIdcG+wHYtILPAkcDqQDiwQkRmqujpgtc3Aiaq6X0TOBibj9FNoVP63NpMPV2Tw2zP60iM5we9wjDHmIF6eESSIyHGVMyIyGvByNhsBbFDVTapaAkwHxgauoKqzVXW/OzsX6OIt7IYjv7iMv7y7ir4dW3KD9R42xkQgr4PXTxWR1u58NnCth+1SgO0B8+kE/7Z/HfBxdW+IyA3ADQCpqQ1r6MZ/fb6eHdmFvDnpWJo1sWfzxpjI42WoykXAYBFJBERVczzuu7ob4dU+fBaRk3ESwXHVva+qk3FuG5GWlnbIB9iRZuWOHKZ+t5krRqZaQTljTMSqMRGISHPgYqA70KTyQaeqBn1GgHMF0DVgvguws5r9DwKmAGer6l5PUTcAZeUV/OntFbRr2Zw/nHWk3+EYY8whebk19B6QgzNkZXEt9r0A6CMiPYAdwDjgisAVRCQVeBu4UlXX12LfEW/anK2s2JHDE1cMpXW8jT1sjIlcXhJBF1U9q7Y7VtUyEbkZ+BSn+ehUVV0lIpPc958G/gq0A55yrzTKVDWttr8r0uzILuSRz9Zxcr/2nHN0J7/DMcaYoLwkgtkicrSq1nowGlX9CPioyrKnA6avB66v7X4jmapy13srUYV7xx5lfQaMMRHPSyI4DrhaRDbj3BoSQK1ncfU+WbmLmWt2c+fP+tO1bQu/wzHGmBp5SQRnhzyKRiK3qJS7ZqxiYOdErhnT3e9wjDHGEy+JoME01/TbQ5+sY09eMVMmptEk1voMGGMaBi+J4EOcZCBAHNADZxzjgSGMq8FZtHU/L83bytWjuzOoS5Lf4RhjjGdeOpQdHTgvIsOAG0MWUQNUWl7B/3t7BUckxvGbM/r5HY4xxtRKre9fqOpi4JgQxNJg/febTazLPMC9Y4+iZXMvF1nGGBM5vPQsviNgNgYYBmSFLKIGZuvefB6b+T1nDTyC0wd09DscY4ypNS9fX1sFTJfhPDN4KzThNCyqyp/fXUnT2BjuPt8emRhjGiYvzwjuARCRVs6s5oU8qgZixrKdfPP9Hu4dO9CGnjTGNFg1PiMQkaNEZAmwElglIotE5KjQhxbZsgtKuPf91QzpmsSEkd38DscYYw6bl1tDk4E7VPVLABE5yV02OnRhRb77P1pLdmEpL110NLExVkbCGNNweRqhrDIJAKjqLLyNUNZozdu0l9cWbuf643vQv1Oi3+EYY0ydeLki2CQif+HHwet/jjPWcFQqLivnT++soGvbeG47ta/f4RhjTJ15uSK4FmiPM27A20AycE0og4pk/5m1kU1Z+fztgqOJbxbrdzjGGFNnQa8IRCQWeENVTwtTPBFtY1YeT325kfMHd+bEvu39DscYY+pF0CsCVS0HCgIGro9qD3y8luZNYvjLuQP8DsUYY+qNl2cERcAKEfkcyK9cqKq3hCyqCLRo634+X53J787sR/tWzf0Oxxhj6o3X6qMfhjqQSKaqPPjJWtq3am7jDBhjGh0vPYunhSOQSDZrfRbzN+/jvguOokUzKypnjGlcvBSdW8FPB6fJARYCf1PVvaEILFJUVCj//GQd3dq1YNwxXf0Oxxhj6p2Xr7cfA+XAK+78OJxBanKA54HzQhJZhHh/+U7WZOTy2LghNLVRx4wxjZCXRDBGVccEzK8Qke9UdYyI/DxUgUWCkrIKHvlsPQM6JXLeoM5+h2OMMSHh5StuSxEZWTkjIiOAlu5sWUiiihDTF2xj274Cfn9WP2KsnpAxppHyckVwPTBVRCpP/geA60UkAbg/ZJH5LL+4jMe/+J5RPdta5zFjTKPmpdXQAuBot1OZqGp2wNuvhyowv039djN78kqYfNWRiNjVgDGm8fIyHkFHEXkWmK6q2SIyQESuC0NsvtmXX8IzX2/izIEdGZbaxu9wjDEmpLw8I3ge+BSofFq6HrgtRPFEhKe+3EBBSRm/O7Of36EYY0zIeUkEyar6OlABoKplOM1JG6Ud2YW8MGcrlwzvQu8OrWrewBhjGjgviSBfRNrhdioTkVE4fQgapX9/vh4EbjvNxhowxkQHL62GfgPMAHqJyHc4YxNcEtKofPJ95gHeWpzOdcf1oHNSvN/hGGNMWHhpNbRIRE4E+uH0KF6nqqUhj8wHD326joRmTfjlSb39DsUYY8LGS6uhZcDvgSJVXdlYk8Cirfv5bHUmN5zQkzYJzfwOxxhjwsbLM4LzcXoQvy4iC0TktyKSGuK4wqqyzHRyy+Zce1wPv8MxxpiwqjERqOpWVf2nqg4HrgAG0cgGr//KLTN9y6m9SWhuZaaNMdHF01lPRLoDlwGX4zQd/X0IYwqrigrlwU/Wkdq2BeOOaVQXOsYY44mX8QjmAU2BN4BLVXVTyKMKo8Ay082aWJlpY0z08XJFMFFV14Y8Eh9Ulpnub2WmjTFRzMtX4AwReVREFrqvR9wCdDUSkbNEZJ2IbBCRP1bzvojI4+77y0VkWK3/BXXwmpWZNsYYT4lgKk7p6cvcVy7wXE0biUgs8CRwNjAAGC8iA6qsdjbQx33dAPzHc+R1lF9cxmNfbGBkj7acZGWmjTFRzMutoV6qenHA/D0istTDdiOADZXPFERkOjAWWB2wzljgBVVVYK6IJIlIJ1XN8Bb+4Xvuu83sySvmmSuHW5lpY0xU83JFUCgix1XOiMgYoNDDdinA9oD5dHdZbddBRG6ovDWVlZXl4VcHtz+/hGe+2sTpAzoyvJuVmTbGRDcvVwQ3AdMqB6YB9gETPWxX3ddsPYx1UNXJwGSAtLS0n7xfW0/N2kC+lZk2xhjAW62hpcBgEUl053M97jsd6Bow3wXYeRjr1Ksd2YVMm7OVi4Z1oW9HKzNtjDFeag21E5HHgVnAlyLymFuWuiYLgD4i0kNEmgHjcKqYBpoBXOW2HhoF5IT6+cBjM9eDwu2nW5lpY4wBb88IpgNZwMU45aezgNdq2sgdwOZmnNHN1gCvq+oqEZkkIpPc1T4CNgEbgP8Cv6z1v6AWvs88wJuL0rny2G6kWJlpY4wBvD0jaKuq9wXM/01ELvCyc1X9COdkH7js6YBpBX7lZV/14eHP1tGiWRN+dbKVmTbGmEpergi+FJFxIhLjvi4DPgx1YPVt8bb9fLrKKTPd1spMG2PMDw55RSAiB3Ba8AhwB/CS+1YMkAfcFfLo6tkJfdtznZWZNsaYgxwyEahqo2pSMyy1DS9cO8LvMIwxJuJ4LUPdBqcMRFzlMlX9OlRBGWOMCR8vZaivB27FaeO/FBgFzAFOCWlkxhhjwsLLw+JbgWOArap6MjAUpwmpMcaYRsBLIihS1SIAEWnujk1gtRmMMaaR8PKMIF1EkoB3gc9FZD8hLgNhjDEmfLzUGrrQnbxbRL4EWgOfhDQqY4wxYeOp1VAlVf0qVIEYY4zxh43WbowxUU6ccj8Nh4hkAVv9jqMGycAev4PwwOKsfw0lVouzfjWEOLuparXj8ja4RNAQiMhCVU3zO46aWJz1r6HEanHWr4YS56HYrSFjjIlylgiMMSbKWSIIjcl+B+CRxVn/GkqsFmf9aihxVsueERhjTJSzKwJjjIlylgiMMSbKWSKoJRHpKiJfisgaEVklIre6y+8WkR0istR9/Sxgmz+JyAYRWSciZ4Yx1i0issKNZ6G7rK2IfC4i37s/20RAnP0CjttSEckVkdsi4ZiKyFQR2S0iKwOW1foYishw97PYICKPi4iEIc6HRGStiCwXkXfcmmGISHcRKQw4rk8HbBPSOIPEWuvP2qdj+lpAjFtEZKm73NdjWmeqaq9avIBOwDB3uhWwHhgA3A38tpr1BwDLgOZAD2AjEBumWLcAyVWW/RP4ozv9R+BBv+OsEl8ssAvoFgnHFDgBGAasrMsxBOYDx+IM/foxcHYY4jwDaOJOPxgQZ/fA9arsJ6RxBom11p+1H8e0yvuPAH+NhGNa15ddEdSSqmao6mJ3+gCwBkgJsslYYLqqFqvqZmAD4OeYmWOBae70NOCCgOWREOepwEZVDdZ7PGyxqjMS375qfr/nYyginYBEVZ2jzpnhhYBtQhanqn6mqmXu7FycwaUOKRxxHirWICLqmFZyv9VfBrwabB/hOqZ1ZYmgDkSkO85APfPcRTe7l+FTA24XpADbAzZLJ3jiqE8KfCYii0TkBndZR1XNACepAR0iIM5A4zj4P1ekHVOo/TFMcaerLg+na3G+jVbqISJLROQrETneXeZ3nLX5rP2O9XggU1W/D1gWicfUE0sEh0lEWgJvAbepai7wH6AXMATIwLlsBOdysKpwtdkdo6rDgLOBX4nICUHW9TNOJwCRZsD5wBvuokg8psEcKi5f4xWRO4Ey4GV3UQaQqqpDgTuAV0QkEX/jrO1n7fffwHgO/sISicfUM0sEh0FEmuIkgZdV9W0AVc1U1XJVrQD+y4+3KtKBrgGbdyFMA/uo6k73527gHTemTPdytfKydbffcQY4G1isqpkQmcfUVdtjmM7Bt2XCFq+ITATOBSa4tyZwb7PsdacX4dx37+tnnIfxWft5TJsAFwGvVS6LxGNaG5YIasm9N/gssEZVHw1Y3ilgtQuBypYGM4BxItJcRHoAfXAeHoU6zgQRaVU5jfPgcKUbz0R3tYnAe37GWcVB37Ii7ZgGqNUxdG8fHRCRUe7fz1UB24SMiJwF/AE4X1ULApa3F5FYd7qnG+cmv+J046jVZ+1nrMBpwFpV/eGWTyQe01rx+2l1Q3sBx+Fc2i0HlrqvnwEvAivc5TOATgHb3InzDWEdYWoxAPTEaW2xDFgF3Okubwd8AXzv/mzrZ5wBv7sFsBdoHbDM92OKk5gygFKcb3fXHc4xBNJwTm4bgSdwe/WHOM4NOPfXK/9On3bXvdj9m1gGLAbOC1ecQWKt9WftxzF1lz8PTKqyrq/HtK4vKzFhjDFRzm4NGWNMlLNEYIwxUc4SgTHGRDlLBMYYE+UsERhjTJSzRGAQkbwa3k8SkV8GzHcWkTdDH1ndicgUERkQgv3eLSK/re/91pWIXC0iT4Rw/51E5AOP64bs70ZEbhaRa+pjX8YSQVQQR10+6yTgh//QqrpTVS+pc2AeVXbUORyqer2qrq7PeOpTXf5tPrkDp+evF0mE7u9mKnBLPe0r6lkiaKTc+uhrROQpnA4uXUXkdyKywC3sdU8127QUkS9EZLFbP32s+9YDQC+3zvpD7r5XutvME5GBAfuYJU799QS3eNgCtxDX2Gp+X4yIPCXOuA4fiMhHInKJ+94WEfmriHwLXCoi492YVorIg+46l4nIo+70rSKyyZ3u5W5XGU+aO50nIn8XkWUiMldEOgasP9eN9d5DXSGJyJ3i1MSfCfQLWN5LRD4Rp7jfNyJyZLD9ishJ4oxp8QqwQkRi3eNa+dncGLDvoJ+Zu841IrJeRL4CxgQsP8/9fJaIyEwR6ege8+9FpH3AZ7BBRJJF5FL3+C4Tka+r+104Hac+OcR7VQX7u7laRN4VkfdFZLP7Df8ON9a5ItI22LFVp6f0FhHxs5Jv4+F3jzZ7heaFUx+9Ahjlzp+BM8C24HwB+AA4wX0vz/3ZBKdkLkAyTs9UoUqt9cB54HbgHne6E7Denf4H8HN3Ogln3IaEKjFeAnzkxnMEsB+4xH1vC/B7d7ozsA1o78b4P5xSvkcAC9x13gQW4FR2nAjc7y6fBaS504rb4xNnTIE/u9MfAOPd6UmVx6NKrMNxer62ABLdY/Nb970vgD7u9Ejgf8H2C5wE5AM93PkbAmJpDizEqb1/yM8sIK5OAcemGfAd8IT7Xht+HJf8euARd/ounGKJuL/jLXd6BZBS+ZlVcwx6AIsC5hOAfwNLcGpZ/Qw4Cnix6t9JNX83V7vHsJUbew5ub13gXwHxVXts3fk7gd/4/X+tMbzsiqBx26qqc93pM9zXEpwrhCNx6qEEEuAfIrIcmIlzUu1Yw+94HbjUnb6MHyuHngH8UZwRnGYBcUBqlW2PA95Q1QpV3QV8WeX9yqJexwCzVDVLnfr6L+OcEHcBLcWpqdQVeAVnMJHjgW+qibUE52QKsAjnxATOoCGVcb9yiH/n8cA7qlqgTrXZGfBDFdrRwBvuv/UZnJNzTfudr059fXCO1VXu9vNwSlj0wdtnNjLg2JQQUAgNp8DZpyKyAvgdUHnlNhWn5g045amfc6e/A54XkV/gDBBUVScgK2B+EPCtOhU3nwR+g1Mq4v1qtq3Ol6p6QFWzcBJB5XYrgO41HFtwiv119vi7TBBN/A7AhFR+wLTgfEt+Jsj6E3C+nQ1X1VIR2YJzAj8kVd0hIntFZBBwOVB5W0OAi1V1XZDNaxqyrzL+YOvNAa7BqUPzDc6J7Vick1JVpep+lQTKqf3ff3X1WGKAbFUdUst9Vf1sfq2qnwauIM6wjDV9ZoeKC+D/gEdVdYaInIQzChiqul1EMkXkFJxEMsFdPklERgLnAEtFZIi6FTVdhQT8PajqnIDpmThfHmqjOGC6ImC+AuezqenYxrkxmTqyK4Lo8SlwrfstCxFJEZEOVdZpDex2k8DJOMNFAhzAuYQ/lOnA73EKxq0I+H2/FnHGZxWRodVs9y1wsXufuiPOLZPqzANOdO9jx+JUKf3Kfe9r4LfuzyXAyUCxquYEibequTj3vsEZGKc6XwMXiki8ewVyHoB7dbBZRC6FHx7MD67FfsE5VjeJU94cEekrTsVYL5/ZPOAkEWnnbn9pwHutgR3u9MQq200BXgJeV9Vyd/+9VHWeqv4V2MPB5Z/Bub3XPci/o6qa/m6CquHYglPmeWW1G5tasUQQJVT1M5zbE3PcWwVv8tP/pC8DaeIMdD8BWOtuuxf4zn2Q+FA1u38T50T3esCy+4CmwHL3AeF91Wz3Fk5Vx5U4l/3zcG4RVI09A/gTzq2jZThjFlSW8v0G54T1tXtC246TYGrjNuAOEZmPc+uhuhgW49x2WerGHXjraQJwnYhUVnqtfDBe435dU4DVwGL3WD2DM9ZwjZ+Ze2zuxrkymolzC6nS3Ti3Vb7BObEHmgG05MfbQgAPiftAHifxLavyu/KBjSLS+xD/joN4+Lvx4lDHFpwH47W9CjHVsOqjxlci0lJV80SkHc6YAmPce//hjKEFUKiqKiLjcB7w/qSVU6Tstz6I05LqX6p6fI0rH7zdhTi3Dv8cmsg8xzEUuENVr/QzjsbCnhEYv30gIkk4LV7uC3cScA0HnnBvY2XjPGeI5P3WiYj8EbgJ99lAbajqO27S9lsy8Be/g2gs7IrAGGOinD0jMMaYKGeJwBhjopwlAmOMiXKWCIwxJspZIjDGmCj3/wFaFTNGDs+YawAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gdd = np.arange(100, 2000, 100)\n",
    "gdd_e = 120  # gdd at emergence\n",
    "gdd_s = 1500  # gdd at spiking\n",
    "gdd_rel = (gdd - gdd_e) / (gdd_s - gdd_e)\n",
    "f_u = 0.12 + 0.88 * np.exp(-2.28 * gdd_rel)\n",
    "f_agb = 1 - f_u\n",
    "plt.plot(gdd, f_agb)\n",
    "\n",
    "plt.xlabel('relative growing degree days (~time)')\n",
    "plt.ylabel('aboveground biomass partitioning fracton')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aade475d-0042-45de-a8b1-35419ab4ffaf",
   "metadata": {},
   "source": [
    "Usually the values of f_agb are taken on annual basis.\n",
    "\n",
    "Reference values for wheat 0.86, potato 0.2, sugarcane 0.8 presented in Table 2-3 in \n",
    "\n",
    "FAO. 2020. WaPOR V2 quality assessment – Technical Report on the Data Quality of the WaPOR FAO Database version 2. Rome.\n",
    "https://doi.org/10.4060/cb2208en"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0de6d6f-25f8-4571-b3ba-973abe038b06",
   "metadata": {},
   "source": [
    "## yield\n",
    "\n",
    "0. accumulate DMP over the growing season\n",
    "1. partition DMP to AGBP (AGB fraction, f_agb) - strictly speaking, a dynamic value but we simplify\n",
    "2. partitioned AGBP to harvestable material (harvest index, HI)\n",
    "3. convert dry yield to fresh yield (moisture content of the harvestable product, $\\theta$, [g DW g$^{-1}$ FW])\n",
    "\n",
    "Crop-specific values of HI and $\\theta$ are from \n",
    "\n",
    "Bastiaanssen, Wim G.M., and Pasquale Steduto. 2017. “The Water Productivity Score (WPS) at Global and Regional Level: Methodology and First Results from Remote Sensing Measurements of Wheat, Rice and Maize.” Science of the Total Environment 575: 595–611. https://doi.org/10.1016/j.scitotenv.2016.09.032."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "26faf45b-b818-452b-98aa-b77510961e16",
   "metadata": {},
   "outputs": [],
   "source": [
    "dmp_season = 100  # kg DM ha-1 season-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "59a45b17-4ceb-400a-83b3-9c6c45347284",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50.588235294117645"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_agb_wheat = 0.86  # from FAO\n",
    "hi_wheat = 0.50 # -\n",
    "theta_wheat = 0.15  # g DW g-1 FW\n",
    "\n",
    "yield_wheat = dmp_season * f_agb_wheat * hi_wheat / (1 - theta_wheat)\n",
    "yield_wheat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "2a7b1267-6e4c-4d06-8ed6-6cca732394b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_agb_maize = 0.90  # guess\n",
    "hi_maize = 0.55 # -\n",
    "theta_maize = 0.26  # g DW g-1 FW\n",
    "\n",
    "yield_maize = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea472ebc-1f2a-4f8c-905f-233cd368b9f7",
   "metadata": {},
   "source": [
    "Find the reference values of f_agb, HI and water content for the remaining study crops: barley, oat, rapeseed\n",
    "\n",
    "Some values are available in \n",
    "- FAO AquaCrop parameterization https://www.fao.org/fileadmin/user_upload/faowater/docs/AquaCropV40Annexes.pdf\n",
    "- WaPOR database methodology https://www.fao.org/3/ca9894en/CA9894EN.pdf\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e52848b-4097-4c52-9295-5408348a2e20",
   "metadata": {},
   "source": [
    "# DMP datasets\n",
    "\n",
    "- Copernicus Global Land Monitoring service https://land.copernicus.eu/global/products/dmp\n",
    "    - global\n",
    "    - 333 m (since 2014), 1000 m (since 1999)\n",
    "    - 10-daily (decadal)\n",
    "    - Sentinel-3 (earlier Proba-V)\n",
    "    - LUE model\n",
    "    - https://notebooks.terrascope.be/ - access to the server for subsetting (otherwise 1 Gb compressed netcdf)\n",
    "- FAO WaPOR - https://wapor.apps.fao.org/catalog/WAPOR_2/1\n",
    "    - Africa + Middle East\n",
    "    - GPP, NPP, ETa\n",
    "    - 250 m - 100 m - 30 m\n",
    "    - 10-daily (decadal)\n",
    "    - MODIS, VIIRS, Landsat, Sentinel-2\n",
    "    - summary statistics on server\n",
    "- MODIS MOD17/MYD17 https://lpdaac.usgs.gov/products/mod17a2hv061/\n",
    "    - global\n",
    "    - 500 m\n",
    "    - 8-daily\n",
    "    - AppEARS subsetting is the most convenient"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "579b8645-14fe-4fe5-bdbb-8eb01da7d51c",
   "metadata": {},
   "source": [
    "# Evapotranspiration (ET)\n",
    "\n",
    "$\\lambda$E - latent heat flux, W m-2\n",
    "ET - evapotranspiration, mm day-1\n",
    "\n",
    "$ET = \\dfrac{\\lambda E}{\\lambda}$,\n",
    "where $\\lambda$ = 2.45 MJ kg H2O\n",
    "\n",
    "1 mm = 1 kg H2O m-2\n",
    "\n",
    "1 W = 1 J s-1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d24a89d4-e8f8-4abb-8ee3-b09f91be546c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lE_W_m2 = 100 # W m-2\n",
    "lambd = 2.45 # MJ kg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d36b4c74-74ba-49a9-9785-cc56a33b5d79",
   "metadata": {},
   "source": [
    "## W to MJ day-1\n",
    "\n",
    "1 W = 1 J s-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0bee8869-bad5-4078-8edc-94205e0092c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.639999999999999"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lE_MJ_m2_s = lE_W_m2 * 1e-6\n",
    "lE_MJ_m2_day = lE_MJ_m2_s * (60*60*24)\n",
    "lE_MJ_m2_day"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe2c7014-4293-4798-8e57-635cc9921e2c",
   "metadata": {},
   "source": [
    "## MJ m-2 to mm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e237a1e0-112d-4ed9-9dd6-68e841d375b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.526530612244897"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ET_mm_day = lE_MJ_m2_day / lambd \n",
    "ET_mm_day"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ef18831-3fc7-4853-9480-f897ff6dcf84",
   "metadata": {},
   "source": [
    "## Shortcut W m-2 -> mm day-1\n",
    "\n",
    "1 W m-2 = 0.035 mm day-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6d92d04e-b6a3-4fe6-a3ae-03d252264bf5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.5000000000000004"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lE_W_m2 * 0.035  # mm day-1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}