leaf_enbalance_eqs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Leaf energy balance equations\n",
    "Based on the following paper:\n",
    "Schymanski, S.J. and Or, D. (2016): [Leaf-scale experiments reveal important omission in the Penman-Monteith equation.](http://www.hydrol-earth-syst-sci-discuss.net/hess-2016-363/) Hydrology and Earth System Sciences Discussions, p.1–33. doi: 10.5194/hess-2016-363.\n",
    "\n",
    "Author: Stan Schymanski (stan.schymanski@env.ethz.ch)\n",
    "\n",
    "Note: This worksheet is prepared for the open source software [sage](http://www.sagemath.org). It relies on definitions provided in Worksheet [Worksheet_setup](Worksheet_setup.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%capture storage\n",
    "# The above redirects all output of the below commands to the variable 'storage' instead of displaying it.\n",
    "# It can be viewed by typing: 'storage()'\n",
    "\n",
    "# Loading modules, settings and custom functions\n",
    "load('temp/Worksheet_setup.sage')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Definition of variables\n",
    "All variables are also listed in a [Table](#Table-of-symbols) at the end of this document."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "v_w"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var2('alpha_a', 'Thermal diffusivity of dry air',meter^2/second)\n",
    "var2('a_s', 'Fraction of one-sided leaf area covered by stomata (1 if stomata are on one side only, 2 if they are on both sides)', meter^2/meter^2)\n",
    "var2('a_sh', 'Fraction of projected area exchanging sensible heat with the air (2)', meter^2/meter^2, value = 2, latexname = 'a_{sh}')\n",
    "var2('c_pa', 'Specific heat of dry air (1010) ', joule/kilogram/kelvin, latexname = 'c_{pa}', value = 1010)\n",
    "#var2('c_pamol', 'Molar specific heat of dry air (29.19) ', joule/mole/kelvin, latexname = 'c_{pa,mol}', value = 29.19) # https://en.wikipedia.org/wiki/Heat_capacity#Specific_heat_capacity\n",
    "var2('C_wa', 'Concentration of water in the free air ',mole/meter^3, latexname = 'C_{wa}')\n",
    "var2('C_wl', 'Concentration of water in the leaf air space ',mole/meter^3, latexname = 'C_{wl}')\n",
    "var2('D_va', 'Binary diffusion coefficient of water vapour in air',meter^2/second, latexname = 'D_{va}')\n",
    "var2('E_lmol', 'Transpiration rate in molar units',mole/second/meter^2,latexname='E_{l,mol}')\n",
    "var2('E_l', 'Latent heat flux from leaf',joule/second/meter^2)\n",
    "var2('epsilon_l', 'Longwave emmissivity of the leaf surface (1.0)', value = 1, units = 1/1)\n",
    "var2('g', 'Gravitational acceleration (9.81)', meter/second^2, value= 9.81)\n",
    "var2('g_bw', 'Boundary layer conductance to water vapour ',meter/second,latexname='g_{bw}')\n",
    "var2('g_bwmol', 'Boundary layer conductance to water vapour ',mole/meter^2/second,latexname='g_{bw,mol}')\n",
    "var2('Gr', 'Grashof number', latexname = 'N_{Gr_L}')\n",
    "var2('g_sw', 'Stomatal conductance to water vapour',meter/second,latexname='g_{sw}')\n",
    "var2('g_swmol', 'Stomatal conductance to water vapour',mole/meter^2/second,latexname='g_{sw,mol}')\n",
    "var2('g_tw', 'Total leaf conductance to water vapour',meter/second,latexname='g_{tw}')\n",
    "var2('g_twmol', 'Total leaf layer conductance to water vapour',mole/meter^2/second,latexname='g_{tw,mol}')\n",
    "var2('h_c', 'Average 1-sided convective transfer coefficient', joule/meter^2/second/kelvin)\n",
    "var2('H_l', 'Sensible heat flux from leaf',joule/second/meter^2)\n",
    "var2('k_a', 'Thermal conductivity of dry air',joule/second/meter/kelvin)\n",
    "var2('lambda_E', 'Latent heat of evaporation (2.45e6)',joule/kilogram,value = 2.45e6)\n",
    "var2('Le', 'Lewis number', latexname = 'N_{Le}')\n",
    "var2('L_l', 'Characteristic length scale for convection (size of leaf)',meter)\n",
    "var2('M_N2', 'Molar mass of nitrogen (0.028)',kilogram/mole,value = 0.028)\n",
    "var2('M_O2', 'Molar mass of oxygen (0.032)',kilogram/mole,value = 0.032)\n",
    "var2('M_w', 'Molar mass of water (0.018)',kilogram/mole,value = 0.018)\n",
    "var2('nu_a', 'Kinematic viscosity of dry air',meter^2/second)\n",
    "var2('Nu', 'Nusselt number',latexname = 'N_{Nu_L}')\n",
    "var2('P_a', 'Air pressure',pascal)\n",
    "var2('Pr', 'Prandtl number (0.71)',latexname = 'N_{Pr}',value = 0.71)\n",
    "var2('P_N2', 'Partial pressure of nitrogen in the atmosphere', pascal, latexname = 'P_{N2}')\n",
    "var2('P_O2', 'Partial pressure of oxygen in the atmosphere', pascal, latexname = 'P_{O2}')\n",
    "var2('P_wa', 'Vapour pressure in the atmosphere', pascal, latexname = 'P_{wa}')\n",
    "var2('P_was', 'Saturation vapour pressure at air temperature', pascal, latexname = 'P_{was}')\n",
    "var2('P_wl', 'Vapour pressure inside the leaf', pascal, latexname = 'P_{wl}')\n",
    "var2('r_bw', 'Boundary layer resistance to water vapour, inverse of $g_{bw}$', second/meter, latexname = 'r_{bw}')\n",
    "var2('r_sw', 'Stomatal resistance to water vapour, inverse of $g_{sw}$', second/meter, latexname = 'r_{sw}')\n",
    "var2('r_tw', 'Total leaf resistance to water vapour, $r_{bv} + r_{sv}$', second/meter, latexname = 'r_{tw}')\n",
    "var2('Re_c', 'Critical Reynolds number for the onset of turbulence',latexname='N_{Re_c}')\n",
    "var2('Re', 'Reynolds number',latexname='N_{Re_L}')\n",
    "var2('rho_a', 'Density of dry air', kilogram/meter^3)\n",
    "var2('rho_al', 'Density of air at the leaf surface', kilogram/meter^3)\n",
    "var2('R_ll', 'Longwave radiation away from leaf',joule/second/meter^2, latexname = 'R_{ll}')\n",
    "var2('R_mol', 'Molar gas constant (8.314472)',joule/mole/kelvin,latexname='R_{mol}', value = 8.314472)\n",
    "var2('R_s', 'Solar shortwave flux',joule/second/meter^2)\n",
    "var2('Sh', 'Sherwood number',latexname = 'N_{Sh_L}')\n",
    "var2('sigm', 'Stefan-Boltzmann constant (5.67e-8)',joule/second/meter^2/kelvin^4,'real',latexname='\\\\sigma',value=5.67e-8)\n",
    "var2('T_a', 'Air temperature',kelvin)\n",
    "var2('T_l', 'Leaf temperature',kelvin)\n",
    "var2('T_w', 'Radiative temperature of objects surrounding the leaf', kelvin)\n",
    "var2('v_w', 'Wind velocity',meter/second)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mathematical derivations\n",
    "## Leaf energy balance\n",
    "The material below follows closely derivations published previously \\citep{schymanski_stomatal_2013}, with re-organisation of equations for greater consistence with the present paper.\n",
    "\n",
    "The leaf energy balance is determined by the dominant energy fluxes between the leaf and its surroundings, including radiative, sensible, and latent energy exchange (linked to mass exchange).\n",
    "\n",
    "\n",
    "In this study we focus on steady-state conditions, in which the energy balance can be written as:\n",
    "##### {eq_Rs_enbal}\n",
    "$$ R_s = R_{ll} + H_l + E_l $$ \n",
    "where $R_s$ absorbed short wave radiation, $R_{ll}$ is the net longwave balance, i.e. the emitted minus the absorbed, $H_l$ is the sensible heat flux away from the leaf and $E_l$ is the latent heat flux away from the leaf. In the above, extensive variables are defined per unit leaf area. Following our previous work \\citep{schymanski_stomatal_2013}, this study considers spatially homogeneous planar leaves, i.e. homogenous illumination and a negligible temperature gradient between the two sides of the leaf. \n",
    "The net longwave emission is represented by the difference between blackbody radiation at leaf temperature ($T_l$) and that at the temperature of the surrounding objects ($T_w$) \\citep{Monteith_principles_2007}:\n",
    "##### {eq_Rll}\n",
    "$$R_{ll} = a_{sH} \\epsilon_l \\sigma (T_{l}^{4} - T_{w}^{4})$$\n",
    "\n",
    "where $a_{sH}$ is the fraction of projected leaf area exchanging radiative and sensible heat (2 for a planar leaf, 1 for a soil surface), $\\epsilon_l$ is the leaf's longwave emmissivity ($\\approx 1$) and $\\sigma$ is the Stefan-Boltzmann constant.\n",
    "Total convective heat transport away from the leaf is represented as:\n",
    "##### {eq_Hl}\n",
    "$$H_l = a_{sH} h_c (T_l - T_a)$$\n",
    "\n",
    "where $h_c$ is the average one-sided convective heat transfer coefficient, determined by properties of the leaf boundary layer.\n",
    "\n",
    "Latent heat flux ($E_l$, W~m$^{-2}$) is directly related to the transpiration rate ($E_{l,mol}$) by:\n",
    "##### {eq_El}\n",
    "$$E_l = E_{l,mol} M_w \\lambda_E$$\n",
    "\n",
    "where $M_w$ is the molar mass of water and $\\lambda_E$ the latent heat of vaporisation. $E_{l,mol}$ (mol~m$^{-2}$~s$^{-1}$) was computed in molar units as a function of the concentration of water vapour within the leaf ($C_{wl}$, mol m$^{-3}$) and in the free air ($C_{wa}$, mol m$^{-3}$) \\citep[Eq. 6.8]{incropera_fundamentals_2006}:\n",
    "##### {eq_Elmol}\n",
    "$$E_{l,mol} = g_{tw}(C_{wl} - C_{wa}) $$\n",
    "\n",
    "where $g_{tw}$ (m~s$^{-1}$) is the total leaf conductance for water vapour, dependent on stomatal ($g_{sw}$) and boundary layer conductance ($g_{bw}$) in the following way:\n",
    "##### {eq_gtw}\n",
    "$$g_{tw} = \\frac{1}{\\frac{1}{g_{sw}} + \\frac{1}{g_{bw}}}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}R_{s} = E_{l} + H_{l} + {R_{ll}}</script></html>"
      ],
      "text/plain": [
       "R_s == E_l + H_l + R_ll"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kilogram/second^3 == kilogram/second^3"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_Rs_enbal = R_s == R_ll + H_l + E_l\n",
    "units_check(eq_Rs_enbal).simplify_full()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{R_{ll}} = {\\left(T_{l}^{4} - T_{w}^{4}\\right)} {a_{sh}} \\epsilon_{l} {\\sigma}</script></html>"
      ],
      "text/plain": [
       "R_ll == (T_l^4 - T_w^4)*a_sh*epsilon_l*sigm"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kilogram/second^3 == kilogram/second^3"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_Rll = R_ll == a_sh*sigm*epsilon_l*(T_l^4 - T_w^4)\n",
    "units_check(eq_Rll).simplify_full()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}H_{l} = -{\\left(T_{a} - T_{l}\\right)} {a_{sh}} h_{c}</script></html>"
      ],
      "text/plain": [
       "H_l == -(T_a - T_l)*a_sh*h_c"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kilogram/second^3 == kilogram/second^3"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_Hl = H_l == a_sh*h_c*(T_l - T_a)\n",
    "units_check(eq_Hl).simplify_full()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}E_{l} = {E_{l,mol}} M_{w} \\lambda_{E}</script></html>"
      ],
      "text/plain": [
       "E_l == E_lmol*M_w*lambda_E"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kilogram/second^3 == kilogram/second^3"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_El = E_l == E_lmol*M_w*lambda_E\n",
    "units_check(eq_El).simplify_full()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}M_{w} \\lambda_{E}</script></html>"
      ],
      "text/plain": [
       "M_w*lambda_E"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44100.0000000000kilogram*meter^2/(mole*second^2)\n"
     ]
    }
   ],
   "source": [
    "ustr = str(units_check(M_w*lambda_E))\n",
    "print str((M_w*lambda_E).subs(cdict)) + ustr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{E_{l,mol}} = -{\\left({C_{wa}} - {C_{wl}}\\right)} {g_{tw}}</script></html>"
      ],
      "text/plain": [
       "E_lmol == -(C_wa - C_wl)*g_tw"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "mole/(meter^2*second) == mole/(meter^2*second)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_Elmol = E_lmol == g_tw*(C_wl - C_wa)\n",
    "units_check(eq_Elmol).simplify_full()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{g_{tw}} = \\frac{{g_{bw}} {g_{sw}}}{{g_{bw}} + {g_{sw}}}</script></html>"
      ],
      "text/plain": [
       "g_tw == g_bw*g_sw/(g_bw + g_sw)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{g_{tw}} = \\frac{1}{\\frac{1}{{g_{bw}}} + \\frac{1}{{g_{sw}}}}</script></html>"
      ],
      "text/plain": [
       "g_tw == (1/(1/g_bw + 1/g_sw))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "meter/second == meter/second"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_gtw = g_tw == 1/(1/g_sw + 1/g_bw)\n",
    "eq_gtw.simplify_full().show()\n",
    "units_check(eq_gtw).simplify_full()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Boundary layer conductance to water vapour\n",
    "The total leaf conductance to water vapour is determined by the boundary layer and stomatal conductances and equal to 1 over the sum of their respectife resistances ($g_{tw} = 1/(r_{sw} + r_{bw})$.\n",
    "The boundary layer conductance for water vapour is equivalent to the mass transfer coefficient for a wet surface \\cite[Eq. 7.41]{incropera_fundamentals_2006}:\n",
    "##### {eq_gbw}\n",
    "$$g_{bw} = N_{Sh_L} D_{va}/L_l$$\n",
    "\n",
    "where $N_{Sh_L}$ is the dimensionless Sherwood number and $D_{va}$ is the diffusivity of water vapour in air.\n",
    "If the convection coefficient for heat is known, the one for mass ($g_{bw}$) can readily be calculated from the relation \\cite[Eq. 6.60]{incropera_fundamentals_2006}:\n",
    "##### {eq_gbw_hc}\n",
    "$$g_{bw} = \\frac{a_s h_c}{\\rho_a c_{pa} N_{Le}^{1-n}}$$\n",
    "\n",
    "where $a_s$ is the fraction of one-sided transpiring surface area in relation to the surface area for sensible heat exchange, $c_{pa}$ is the constant-pressure heat capacity of air, $n$ is an empirical constant ($n=1/3$ for general purposes) and $N_{Le}$ is the dimensionless Lewis number, defined as \\cite[Eq. 6.57]{incropera_fundamentals_2006}:\n",
    "##### {eq_Le}\n",
    "$$N_{Le} = \\alpha_a/D_{va}$$\n",
    "\n",
    "where $\\alpha_a$ is the thermal diffusivity of air. The value of $a_s$ was set to 1 for leaves with stomata on one side only, and to 2 for stomata on both sides. Other values could be used for leaves only partly covered by stomata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{g_{bw}} = \\frac{{D_{va}} {N_{Sh_L}}}{L_{l}}</script></html>"
      ],
      "text/plain": [
       "g_bw == D_va*Sh/L_l"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{g_{bw}} = \\frac{a_{s} h_{c}}{{N_{Le}}^{\\frac{2}{3}} {c_{pa}} \\rho_{a}}</script></html>"
      ],
      "text/plain": [
       "g_bw == a_s*h_c/(Le^(2/3)*c_pa*rho_a)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "meter/second == meter/second\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{N_{Le}} = \\frac{\\alpha_{a}}{{D_{va}}}</script></html>"
      ],
      "text/plain": [
       "Le == alpha_a/D_va"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1 == 1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_gbw = g_bw == D_va*Sh/L_l\n",
    "units_check(eq_gbw).simplify_full()\n",
    "eq_gbw_hc = g_bw == a_s*h_c/(rho_a*c_pa*Le^(2/3))\n",
    "print units_check(eq_gbw_hc).simplify_full()\n",
    "eq_Le = Le == alpha_a/D_va\n",
    "units_check(eq_Le).simplify_full()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Effect of leaf temperature on the leaf-air vapour concentration gradient\n",
    "\\label{sec_Tleaf}\n",
    "The concentration difference in Eq. [eq_Elmol](#{eq_Elmol}) is a function of the temperature and the vapour pressure differences between the leaf and the free air. Assuming that water vapour behaves like an ideal gas, we can express its concentration as:\n",
    "##### {eq_Cwl}\n",
    "$$C_{wl} = \\frac{P_{wl}}{R_{mol} T_l}\n",
    "$$\n",
    "\n",
    "where $P_{wl}$ is the vapour pressure inside the leaf, $R_{mol}$ is the universal gas constant and $T_l$ is leaf temperature. A similar relation holds for the vapour concentration in free air, $C_{wa} = P_{wa}/(R_{mol} T_l)$. In this study the vapour pressure inside the leaf is assumed to be the saturation vapour pressure at leaf temperature, which is computed using the Clausius-Clapeyron relation \\citep[Eq. B.3]{hartmann_global_1994}:\n",
    "##### {eq_Pwl}\n",
    "$$P_{wl} = 611 \\exp \\left( \\frac{\\lambda_E M_w}{R_{mol}} \\left( \\frac{1}{273} - \\frac{1}{T_l} \\right) \\right)\n",
    "$$\n",
    "where $\\lambda_E$ is the latent heat of vaporisation and $M_w$ is the molar mass of water. \n",
    "\n",
    "Note that the dependence of the leaf-air water concentration difference ($C_{wl} - C_{wa}$) in Eq. [eq_Cwl](#{eq_Cwl}) is very sensitive to leaf temperature. For example if the leaf temperature increases by 5K relative to air temperature, $C_{wl} - C_{wa}$ would double, while if leaf temperature decreased by 6K, $C_{wl} - C_{wa}$ would go to 0 at 70\\% relative humidity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p> </p>\n",
    "<p> </p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{C_{wl}} = \\frac{{P_{wl}}}{{R_{mol}} T_{l}}</script></html>"
      ],
      "text/plain": [
       "C_wl == P_wl/(R_mol*T_l)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mole/meter^3 == mole/meter^3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{P_{wl}} = 611 \\, e^{\\left(-\\frac{M_{w} \\lambda_{E} {\\left(\\frac{273}{T_{l}} - 1\\right)}}{273 \\, {R_{mol}}}\\right)}</script></html>"
      ],
      "text/plain": [
       "P_wl == 611*e^(-1/273*M_w*lambda_E*(273/T_l - 1)/R_mol)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eq_Cwl = C_wl == P_wl/(R_mol*T_l)\n",
    "print units_check(eq_Cwl).simplify_full()\n",
    "eq_Pwl = P_wl == 611*exp(lambda_E*M_w/R_mol*(1/273 - 1/T_l))\n",
    "show(eq_Pwl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Concentration or vapour pressure gradient driving transpiration?\n",
    "\n",
    "Note that $E_{l,mol}$ is commonly expressed as a function of the vapour pressure difference between the free air ($P_{wa}$) and the leaf ($P_{wl}$), in which the conductance ($g_{tw,mol}$) is expressed in molar units (mol~m$^{-2}$~s$^{-1}$):\n",
    "##### {eq_Elmol_conv}\n",
    "$$ E_{l,mol} =  g_{tw,mol} \\frac{P_{wl} - P_{wa}}{P_a}$$\n",
    "\n",
    "For $P_{wl} = P_{wa}$, Eq. [eq_Elmol](#{eq_Elmol}) can still give a flux, whereas Eq. [eq_Elmol_conv](#{eq_Elmol_conv}) gives zero flux. This is because the concentrations of vapour in air (mol~m$^{-3}$) can differ due to differences in tempertaure, even if the partial vapour pressures are the same (see Eq. [eq_Cwl](#{eq_Cwl})). Therefore, the relation between $g_{tw}$ and $g_{v,mol}$ has an asymptote at the equivalent temperature. It can be obtained by combining Eqs. [eq_Elmol](#{eq_Elmol}) and [eq_Elmol_conv](#{eq_Elmol_conv}) and solving for $g_{tw,mol}$:\n",
    "##### {eq_gtwmol_gtw}\n",
    "$$g_{tw,mol} = g_{tw}\\frac{P_a(P_{wa} T_l - P_{wl} T_a)}{(P_{wa} - P_{wl}) R_{mol} T_a T_l)}$$\n",
    "\n",
    "For $T_l = T_a$, the relation simplifies to:\n",
    "##### {eq_gtwmol_gtw_iso}\n",
    "$$g_{tw,mol} = g_{tw}\\frac{P_a }{R_{mol} T_a}$$\n",
    "\n",
    "which, for typical values of $P_a$ and $T_a$ amounts to $g_{tw,mol}\\approx40$~mol~m$^{-3}g_{tw}$.  For all practical purposes, we found that Eqs. [eq_Elmol](#{eq_Elmol}) and [eq_Elmol_conv](#{eq_Elmol_conv}) with $g_{tw,mol} = g_{tw}\\frac{P_a }{R_{mol} T_a}$ give similar results when plotted as functions of leaf temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{E_{l,mol}} = -\\frac{{\\left({P_{wa}} - {P_{wl}}\\right)} {g_{tw,mol}}}{P_{a}}</script></html>"
      ],
      "text/plain": [
       "E_lmol == -(P_wa - P_wl)*g_twmol/P_a"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mole/(meter^2*second) == mole/(meter^2*second)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{g_{tw,mol}} = -\\frac{{\\left(P_{a} {P_{wl}} T_{a} - P_{a} {P_{wa}} T_{l}\\right)} {g_{tw}}}{{\\left({P_{wa}} - {P_{wl}}\\right)} {R_{mol}} T_{a} T_{l}}</script></html>"
      ],
      "text/plain": [
       "g_twmol == -(P_a*P_wl*T_a - P_a*P_wa*T_l)*g_tw/((P_wa - P_wl)*R_mol*T_a*T_l)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mole/(meter^2*second) == mole/(meter^2*second)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{g_{tw,mol}} = \\frac{P_{a} {g_{tw}}}{{R_{mol}} T_{a}}</script></html>"
      ],
      "text/plain": [
       "g_twmol == P_a*g_tw/(R_mol*T_a)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mole/(meter^2*second) == mole/(meter^2*second)\n"
     ]
    }
   ],
   "source": [
    "eq_Elmol_conv = E_lmol == g_twmol*(P_wl-P_wa)/P_a\n",
    "print units_check(eq_Elmol_conv).simplify_full()\n",
    "soln = solve([eq_Elmol.subs(eq_Cwl, eq_Cwl(C_wl = C_wa, P_wl = P_wa, T_l = T_a)), eq_Elmol_conv], g_twmol, E_lmol)\n",
    "eq_gtwmol_gtw = soln[0][0]\n",
    "print units_check(eq_gtwmol_gtw).simplify_full()\n",
    "eq_gtwmol_gtw_iso = eq_gtwmol_gtw(T_l = T_a).simplify_full()\n",
    "print units_check(eq_gtwmol_gtw_iso).simplify_full()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Model closure\n",
    "\n",
    "Given climatic forcing as $P_a$, $T_a$, $R_s$, $P_{wa}$ and $v_w$, and leaf-specific parameters $a_s$, $a_{sH}$, $L_l$ and $g_{sw}$, we need to compute $C_{wa}$, $h_c$, $g_{bw}$ and a series of other derived variables, as described below. \n",
    "\n",
    "The vapour concentration in the free air can be computed from vapour pressure analogously to Eq. [eq_Cwl](#{eq_Cwl}):\n",
    "##### {eq_Cwa}\n",
    "$$C_{wa} = \\frac{P_{wa}}{R_{mol} T_a}$$\n",
    "\n",
    "The heat transfer coefficient ($h_c$) for a flat plate can be determined using the non-dimensional Nusselt number ($N_{Nu_L}$):\n",
    "##### {eq_hc}\n",
    "$$h_{c} = k_a\\frac{N_{Nu_L}}{L_l}$$\n",
    "\n",
    "where $k_a$ is the thermal conductivity of the air in the boundary layer and $L_l$ is a characteristic length scale of the leaf. \n",
    "\n",
    "For sufficiently high wind speeds, inertial forces drive the convective heat transport (forced convection) and the relevant dimensionless number is the Reynolds number ($N_{Re_L}$), which defines the balance between inertial and viscous forces \\cite[Eq. 6.41]{incropera_fundamentals_2006}:\n",
    "##### {eq_Re}\n",
    "$$N_{Re_L} = \\frac{v_w L_l}{\\nu_a}$$\n",
    "\n",
    "where $v_w$ is the wind velocity (m s$^{-1}$), $\\nu_a$ is the kinematic viscosity of the air and $L_l$ is taken as the length of the leaf in wind direction.\n",
    "\n",
    "In the absence of wind, buoyancy forces, driven by the density gradient between the air at the surface of the leaf and the free air dominate convective heat exchange (``free'' or ``natural convection''). The relevant dimensionless number here is the Grashof number ($N_{Gr_L}$), which defines the balance between buoyancy and viscous forces \\cite[Eqs. 9.3 and 9.65]{incropera_fundamentals_2006}:\n",
    "##### {eq_Gr}\n",
    "$$N_{Gr_L} =\\frac{g(\\frac{\\rho_a - \\rho_{al}}{\\rho_{al}}) L_l^3}{\\nu_a^2}$$\n",
    "\n",
    "where $g$ is gravity, while $\\rho_a$ and $\\rho_{al}$ are the densities of the gas in the atmosphere and at the leaf surface respectively. \n",
    "\n",
    "For $N_{Gr_L} \\ll N_{Re_L}^2$, forced convection is dominant and free convection can be neglected, whereas for $N_{Gr_L} \\gg N_{Re_L}^2$ free convection is dominant and forced convection can be neglected \\cite[P. 565]{incropera_fundamentals_2006}. For simplicity, the analysis is limited to forced conditions, which is satisfied by considering wind speeds greater than 0.5~m~s$^{-1}$ for $5 \\times 5$ cm leaves.  %See leaf_chamber4.sws or leaf_capacitance_steady_state3.sws.\n",
    "\n",
    "The average Nusselt number under forced convection was calculated as a function of the average Reynolds number and a critical Reynolds number ($N_{Re_c}$) that determines the onset of turbulence and depends on the level of turbulence in the free air stream or leaf surface properties \\cite[P. 412]{incropera_fundamentals_2006}\n",
    "\n",
    "##### {eq_Nu_forced_all}\n",
    "$$N_{Nu_L} = (0.037 N_{Re_L}^{4/5} - C_1)N_{Pr}^{1/3}$$\n",
    "\n",
    "with\n",
    "##### {eq_C1}\n",
    "$$C_1 = 0.037 C_2^{4/5} - 0.664 C_2^{1/2}$$\n",
    "\n",
    "and\n",
    "##### {eq_C2}\n",
    "$$C_2 = \\frac{N_{Re_L} + N_{Re_c} - |N_{Re_c} - N_{Re_L}|}{2}$$\n",
    "\n",
    "Eq. [eq_C2](#{eq_C2}) was introduced to make Eq.  [eq_Nu_forced_all](#{eq_Nu_forced_all}) valid for all Reynolds numbers, and following considerations explained in our previous work \\citep{schymanski_stomatal_2013}, we chose $N_{Re_c}=3000$ in the present simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{C_{wa}} = \\frac{{P_{wa}}}{{R_{mol}} T_{a}}</script></html>"
      ],
      "text/plain": [
       "C_wa == P_wa/(R_mol*T_a)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mole/meter^3 == mole/meter^3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}h_{c} = \\frac{{N_{Nu_L}} k_{a}}{L_{l}}</script></html>"
      ],
      "text/plain": [
       "h_c == Nu*k_a/L_l"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kilogram/(kelvin*second^3) == kilogram/(kelvin*second^3)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{N_{Re_L}} = \\frac{L_{l} v_{w}}{\\nu_{a}}</script></html>"
      ],
      "text/plain": [
       "Re == L_l*v_w/nu_a"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 == 1\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{N_{Gr_L}} = \\frac{L_{l}^{3} g {\\left(\\rho_{a} - \\rho_{\\mathit{al}}\\right)}}{\\nu_{a}^{2} \\rho_{\\mathit{al}}}</script></html>"
      ],
      "text/plain": [
       "Gr == L_l^3*g*(rho_a - rho_al)/(nu_a^2*rho_al)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 == 1\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{N_{Nu_L}} = {N_{Pr}}^{\\frac{1}{3}} {\\left(-0.0370000000000000 \\, {\\left({N_{Re_L}} + {N_{Re_c}} - \\frac{1}{2} \\, {\\left| {N_{Re_L}} - {N_{Re_c}} \\right|}\\right)}^{\\frac{4}{5}} + 0.0370000000000000 \\, {N_{Re_L}}^{\\frac{4}{5}} + 0.664000000000000 \\, \\sqrt{{N_{Re_L}} + {N_{Re_c}} - \\frac{1}{2} \\, {\\left| {N_{Re_L}} - {N_{Re_c}} \\right|}}\\right)}</script></html>"
      ],
      "text/plain": [
       "Nu == Pr^(1/3)*(-0.0370000000000000*(Re + Re_c - 1/2*abs(Re - Re_c))^(4/5) + 0.0370000000000000*Re^(4/5) + 0.664000000000000*sqrt(Re + Re_c - 1/2*abs(Re - Re_c)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eq_Cwa = C_wa == P_wa/(R_mol*T_a)\n",
    "print units_check(eq_Cwa)\n",
    "eq_hc = h_c == k_a*Nu/L_l\n",
    "print units_check(eq_hc)\n",
    "eq_Re = Re == v_w*L_l/nu_a\n",
    "print units_check(eq_Re)\n",
    "eq_Gr = Gr == g*(rho_a - rho_al)/rho_al * L_l^3/nu_a^2\n",
    "print units_check(eq_Gr)\n",
    "C2 = (Re + Re_c - (abs(Re - Re_c))/2)\n",
    "C1 = 0.037*C2^(4/5) - 0.664*C2^(1/2)\n",
    "eq_Nu_forced_all = Nu == (0.037*Re^(4/5) - C1)*Pr^(1/3)\n",
    "eq_Nu_forced_all.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### Thermodynamic variables\n",
    "In order to simulate steady state leaf temperatures and the leaf energy balance terms using the above equations, it is necessary to calculate $\\rho_a$, $D_{wa}$, $\\alpha_a$, $k_a$, and $\\nu_a$, while $L_l$, $Re_c$ and $g_{sv}$ are input parameters, and $P_{wa}$ and $v_w$ (vapour pressure and wind speed) are part of the environmental forcing.\n",
    "$D_{wa}$, $\\alpha_a$, $k_a$ and $\\nu_a$ were parameterised as functions of air temperature ($T_a$) only, by fitting linear curves to published data \\cite[Table A.3]{Monteith_principles_2007}:\n",
    "##### {eq_Dwa}\n",
    "$$D_{wa} = (1.49\\times 10^{-7})T_a - 1.96\\times 10^{-5}$$\n",
    "\n",
    "##### {eq_alphaa}\n",
    "$$\\alpha_a = (1.32\\times 10^{-7})T_a - 1.73\\times 10^{-5}$$\n",
    "\n",
    "##### {eq_ka}\n",
    "$$k_a = (6.84\\times 10^{-5})T_a + 5.62\\times 10^{-3}$$\n",
    " \n",
    "##### {eq_nua}\n",
    "$$\\nu_a = (9\\times 10^{-8})T_a - 1.13\\times 10^{-5}$$\n",
    "\n",
    "\n",
    "Assuming that air and water vapour behave like an ideal gas, and that dry air is composed of 79\\% N$_2$ and 21\\% O$_2$, we calculated air density as a function of temperature, vapour pressure and the partial pressures of the other two components using the ideal gas law:\n",
    "##### {eq_rhoa_Pa_Ta}\n",
    "$$\\rho_a = \\frac{n_a M_a}{V_a} = M_a \\frac{P_a}{R_{mol} T_a}$$\n",
    "\n",
    "where $n_a$ is the amount of matter (mol), $M_a$ is the molar mass (kg~mol$^{-1}$), $P_a$ the pressure, $T_a$ the temperature and $R_{mol}$ the molar universal gas constant.\n",
    "This equation was used for each component, i.e. water vapour, N$_2$ and O$_2$, where the partial pressures of N$_2$ and O$_2$ are calculated from atmospheric pressure minus vapour pressure, yielding:\n",
    "##### {eq_rhoa_Pwa_Ta}\n",
    "$$\\rho_a = \\frac{M_w P_{wa} + M_{N_2} P_{N_2} + M_{O_2} P_{O_2}}{R_{mol} T_a}$$\n",
    "\n",
    "where $M_{N_2}$ and $M_{O_2}$ are the molar masses of nitrogen and oxygen respectively, while $P_{N_2}$ and $P_{O_2}$ are their partial pressures, calculated as:\n",
    "\n",
    "##### {eq_PN2}\n",
    "$$P_{N_2} = 0.79(P_a - P_{wa})$$\n",
    "\n",
    "and\n",
    "##### {eq_PO2}\n",
    "$$P_{O_2} = 0.21(P_a - P_{wa})$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\rho_{a} = \\frac{M_{N_{2}} {P_{N2}} + M_{O_{2}} {P_{O2}} + M_{w} {P_{wa}}}{{R_{mol}} T_{a}}</script></html>"
      ],
      "text/plain": [
       "rho_a == (M_N2*P_N2 + M_O2*P_O2 + M_w*P_wa)/(R_mol*T_a)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kilogram/meter^3 == kilogram/meter^3\n"
     ]
    }
   ],
   "source": [
    "eq_Dva = D_va == 1.49e-07*T_a - 1.96e-05\n",
    "eq_alphaa = alpha_a == (1.32e-07)*T_a - 1.73e-05\n",
    "eq_ka = k_a == 6.84e-05*T_a + 5.63e-03\n",
    "eq_nua = nu_a == 9.e-08*T_a - 1.13e-05\n",
    "eq_rhoa_Pwa_Ta = rho_a == (M_w*P_wa + M_N2*P_N2 + M_O2*P_O2)/(R_mol*T_a)\n",
    "print units_check(eq_rhoa_Pwa_Ta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}P_{a} = {P_{N2}} + {P_{O2}} + {P_{wa}}</script></html>"
      ],
      "text/plain": [
       "P_a == P_N2 + P_O2 + P_wa"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kilogram/(meter*second^2) == kilogram/(meter*second^2)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{P_{O2}} = \\frac{21}{100} \\, P_{a} - \\frac{21}{100} \\, {P_{wa}}</script></html>"
      ],
      "text/plain": [
       "P_O2 == 21/100*P_a - 21/100*P_wa"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kilogram/(meter*second^2) == kilogram/(meter*second^2)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}{P_{N2}} = \\frac{79}{100} \\, P_{a} - \\frac{79}{100} \\, {P_{wa}}</script></html>"
      ],
      "text/plain": [
       "P_N2 == 79/100*P_a - 79/100*P_wa"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kilogram/(meter*second^2) == kilogram/(meter*second^2)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_Pa = P_a == P_N2 + P_O2 + P_wa\n",
    "print units_check(eq_Pa)\n",
    "eq_PN2_PO2 = P_N2 == 0.79/0.21*P_O2\n",
    "soln = solve([eq_Pa, eq_PN2_PO2], P_O2, P_N2)\n",
    "eq_PO2 = soln[0][0]\n",
    "print units_check(eq_PO2)\n",
    "eq_PN2 = soln[0][1]\n",
    "units_check(eq_PN2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\rho_{a} = \\frac{79 \\, M_{N_{2}} {\\left(P_{a} - {P_{wa}}\\right)} + 21 \\, M_{O_{2}} {\\left(P_{a} - {P_{wa}}\\right)} + 100 \\, M_{w} {P_{wa}}}{100 \\, {R_{mol}} T_{a}}</script></html>"
      ],
      "text/plain": [
       "rho_a == 1/100*(79*M_N2*(P_a - P_wa) + 21*M_O2*(P_a - P_wa) + 100*M_w*P_wa)/(R_mol*T_a)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kilogram/meter^3 == kilogram/meter^3"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_rhoa = eq_rhoa_Pwa_Ta.subs(eq_PN2, eq_PO2)\n",
    "units_check(eq_rhoa).simplify_full()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(C_wa, R_s, T_a, T_l, T_w, g_tw, h_c)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Energy balance as a function of variables that are independent of leaf temperature\n",
    "eqenbalTl = (eq_Rs_enbal - R_s).rhs().subs(eq_El, eq_Hl, eq_Rll).subs(eq_Elmol).subs(eq_Cwl).subs(eq_Pwl)\n",
    "eqenbalTl.subs(cdict).args()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Below, we create a copy of cdict, which is a dictionary with general constants, generated during the definition of variables at the top of the worksheet, and some example values for external forcing and leaf properties needed to compute steady-state fluxes and leaf temperature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Input values\n",
    "vdict = cdict.copy()\n",
    "vdict[a_s] = 1.0    # one sided stomata\n",
    "vdict[g_sw] = 0.01    \n",
    "vdict[T_a] = 273 + 25.5\n",
    "vdict[T_w] = vdict[T_a] # Wall temperature equal to air temperature\n",
    "vdict[P_a] = 101325\n",
    "rha = 1\n",
    "vdict[P_wa] = rha*eq_Pwl.rhs()(T_l = T_a).subs(vdict)\n",
    "vdict[L_l] = 0.03\n",
    "vdict[Re_c] = 3000\n",
    "vdict[R_s] = 600\n",
    "vdict[v_w] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we compute all derived variables, using the equations described above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "26.1863624980041"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Nusselt number\n",
    "vdict[nu_a] = eq_nua.rhs().subs(vdict)\n",
    "vdict[Re] = eq_Re.rhs().subs(vdict)\n",
    "vdict[Nu] = eq_Nu_forced_all.rhs().subs(vdict)\n",
    "vdict[Nu]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22.7362219510171"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# h_c\n",
    "vdict[k_a] = eq_ka.rhs().subs(vdict)\n",
    "vdict[h_c] = eq_hc.rhs().subs(vdict)\n",
    "vdict[h_c]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0209367439791524"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# gbw\n",
    "vdict[D_va] = eq_Dva.rhs().subs(vdict)\n",
    "vdict[alpha_a] = eq_alphaa.rhs().subs(vdict)\n",
    "vdict[rho_a] =  eq_rhoa.rhs().subs(vdict)\n",
    "vdict[Le] =  eq_Le.rhs().subs(vdict)\n",
    "vdict[g_bw] = eq_gbw_hc.rhs().subs(vdict)\n",
    "vdict[g_bw]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1.13400000000000e-7)*T_l^4 - 900.306522304087\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "45.4724439020342*T_l - 13573.5245047572"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Hl, Rll\n",
    "vdict[R_ll] = eq_Rll.rhs().subs(vdict)\n",
    "print vdict[R_ll]\n",
    "vdict[H_l] = eq_Hl.rhs().subs(vdict)\n",
    "vdict[H_l]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21932.0720390026*e^(-5304.00487246815/T_l + 19.4285892764401)/T_l - 386.319257853523"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# El\n",
    "vdict[g_tw] =  eq_gtw.rhs().subs(vdict)\n",
    "vdict[C_wa] = eq_Cwl.rhs()(P_wl = P_wa, T_l = T_a).subs(vdict)\n",
    "vdict[P_wl] = eq_Pwl.rhs().subs(vdict)\n",
    "vdict[C_wl] = eq_Cwl.rhs().subs(vdict)\n",
    "vdict[E_lmol] = eq_Elmol.rhs().subs(vdict)\n",
    "vdict[E_l] = eq_El.rhs().subs(eq_Elmol).subs(vdict)\n",
    "vdict[E_l]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have a dictionary with $E_l$, $H_l$ and $R_{ll}$ as functions of leaf temperature ($T_l$) only. At steady state, [eq_Rs_enbal](#{eq_Rs_enbal}) must be satisfied, so we will substitute the above dictionary into eq_Rs_enbal, subtract $R_s$ and do a numerical search for the root of the equation to obtain steady-state leaf temperature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "305.6506484227355"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Numerical search for steady-state Tl\n",
    "vdict[T_l] = find_root((eq_Rs_enbal - R_s).rhs().subs(vdict), 273, 373)\n",
    "vdict[T_l]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we define a function that performs all the above calculations given a dictionary with the instantaneous forcing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fun_SS(vdict1):\n",
    "    '''\n",
    "    Steady-state T_l, R_ll, H_l and E_l under forced conditions.\n",
    "    Parameters are given in a dictionary (vdict) with the following entries:\n",
    "    a_s, a_sh, L_l, P_a, P_wa, R_s, Re_c, T_a, g_sw, v_w\n",
    "    ''' \n",
    "    vdict = vdict1.copy()\n",
    "\n",
    "    # Nusselt number\n",
    "    vdict[nu_a] = eq_nua.rhs().subs(vdict)\n",
    "    vdict[Re] = eq_Re.rhs().subs(vdict)\n",
    "    vdict[Nu] = eq_Nu_forced_all.rhs().subs(vdict)\n",
    "    \n",
    "    # h_c\n",
    "    vdict[k_a] = eq_ka.rhs().subs(vdict)\n",
    "    vdict[h_c] = eq_hc.rhs().subs(vdict)\n",
    " \n",
    "    # gbw\n",
    "    vdict[D_va] = eq_Dva.rhs().subs(vdict)\n",
    "    vdict[alpha_a] = eq_alphaa.rhs().subs(vdict)\n",
    "    vdict[rho_a] =  eq_rhoa.rhs().subs(vdict)\n",
    "    vdict[Le] =  eq_Le.rhs().subs(vdict)\n",
    "    vdict[g_bw] = eq_gbw_hc.rhs().subs(vdict)   \n",
    "    \n",
    "    # Hl, Rll\n",
    "    vdict[R_ll] = eq_Rll.rhs().subs(vdict)\n",
    "    vdict[H_l] = eq_Hl.rhs().subs(vdict)   \n",
    "\n",
    "    # El\n",
    "    vdict[g_tw] =  eq_gtw.rhs().subs(vdict)\n",
    "    vdict[C_wa] = eq_Cwl.rhs()(P_wl = P_wa, T_l = T_a).subs(vdict)\n",
    "    vdict[P_wl] = eq_Pwl.rhs().subs(vdict)\n",
    "    vdict[C_wl] = eq_Cwl.rhs().subs(vdict)\n",
    "    vdict[E_lmol] = eq_Elmol.rhs().subs(vdict)\n",
    "    vdict[E_l] = eq_El.rhs().subs(eq_Elmol).subs(vdict)    \n",
    "\n",
    "    # Tl\n",
    "    try:\n",
    "        vdict[T_l] = find_root((eq_Rs_enbal - R_s).rhs().subs(vdict), 273, 373)\n",
    "    except: \n",
    "        print 'too many unknowns for finding T_l: ' + str((eq_Rs_enbal - R_s).rhs().subs(vdict).args())\n",
    "    \n",
    "    # Re-inserting T_l\n",
    "    Tlss = vdict[T_l]\n",
    "    for name1 in [C_wl, P_wl, R_ll, H_l, E_l, E_lmol]:\n",
    "        vdict[name1] = vdict[name1].subs(T_l = Tlss)\n",
    "    \n",
    "    # Test for steady state\n",
    "    if n((E_l + H_l + R_ll - R_s).subs(vdict))>1.:\n",
    "        return 'error in energy balance: El + Hl + Rll - R_s = ' + str(n((E_l + H_l + R_ll - R_s).subs(vdict))) \n",
    "    return vdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C_wa            1.29441408346663\n",
      "C_wl            1.91570361006325\n",
      "D_va            0.0000248765000000000\n",
      "E_l             185.424519010311\n",
      "E_lmol          0.00420463761928142\n",
      "H_l             325.157459266011\n",
      "L_l             0.0300000000000000\n",
      "Le              0.888469037042992\n",
      "M_N2            0.0280000000000000\n",
      "M_O2            0.0320000000000000\n",
      "M_w             0.0180000000000000\n",
      "Nu              26.1863624980041\n",
      "P_a             101325         \n",
      "P_wa            3212.56734153661\n",
      "P_wl            4868.42309771766\n",
      "Pr              0.710000000000000\n",
      "R_ll            89.4180217236781\n",
      "R_mol           8.31447200000000\n",
      "R_s             600            \n",
      "Re              1927.40122068744\n",
      "Re_c            3000           \n",
      "T_a             298.500000000000\n",
      "T_l             305.650648423  \n",
      "T_w             298.500000000000\n",
      "a_s             1.00000000000000\n",
      "a_sh            2              \n",
      "alpha_a         0.0000221020000000000\n",
      "c_pa            1010           \n",
      "epsilon_l       1              \n",
      "g               9.81000000000000\n",
      "g_bw            0.0209367439791524\n",
      "g_sw            0.0100000000000000\n",
      "g_tw            0.00676759777734245\n",
      "h_c             22.7362219510171\n",
      "k_a             0.0260474000000000\n",
      "lambda_E        2.45000000000000e6\n",
      "nu_a            0.0000155650000000000\n",
      "rho_a           1.16339248053449\n",
      "sigm            5.67000000000000e-8\n",
      "v_w             1              \n"
     ]
    }
   ],
   "source": [
    "# Test\n",
    "vdict = cdict.copy()\n",
    "#vdict= {}\n",
    "vdict[a_s] = 1.0    # one sided stomata\n",
    "vdict[g_sw] = 0.01    \n",
    "vdict[T_a] = 273 + 25.5\n",
    "vdict[T_w] = vdict[T_a] # Wall temperature equal to air temperature\n",
    "vdict[P_a] = 101325\n",
    "rha = 1\n",
    "vdict[P_wa] = rha*eq_Pwl.rhs()(T_l = T_a).subs(vdict)\n",
    "vdict[L_l] = 0.03\n",
    "vdict[Re_c] = 3000\n",
    "vdict[R_s] = 600\n",
    "vdict[v_w] = 1\n",
    "resdict = fun_SS(vdict)\n",
    "\n",
    "fun_dict_print(resdict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### Calculation using known $T_l$\n",
    "If leaf temperature ($T_l$) is known (e.g. measured), but stomatal conductance $g_{sw}$ is not known, we can still calculate $H_l$ and $R_{ll}$, and then obtain $E_l$ from the energy balance equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}E_{l} = -H_{l} - {R_{ll}} + R_{s}</script></html>"
      ],
      "text/plain": [
       "E_l == -H_l - R_ll + R_s"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kilogram/second^3 == kilogram/second^3"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_El_enbal = solve(eq_Rs_enbal, E_l)[0]\n",
    "units_check(eq_El_enbal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "E_l == 185.462402956757"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test with known Tl but no g_sw\n",
    "vdict = cdict.copy()\n",
    "vdict[a_s] = 1.0    # one sided stomata\n",
    "   \n",
    "vdict[T_a] = 273 + 25.5\n",
    "vdict[T_w] = vdict[T_a] # Wall temperature equal to air temperature\n",
    "vdict[P_a] = 101325\n",
    "rha = 1\n",
    "vdict[P_wa] = rha*eq_Pwl.rhs()(T_l = T_a).subs(vdict)\n",
    "vdict[L_l] = 0.03\n",
    "vdict[Re_c] = 3000\n",
    "vdict[R_s] = 600\n",
    "vdict[v_w] = 1\n",
    "\n",
    "vdict[C_wa] = eq_Cwl.rhs()(P_wl = P_wa, T_l = T_a).subs(vdict)\n",
    "vdict[nu_a] = eq_nua.rhs().subs(vdict)\n",
    "vdict[Re] = eq_Re.rhs().subs(vdict)\n",
    "vdict[Nu] = eq_Nu_forced_all.rhs().subs(vdict)\n",
    "vdict[k_a] = eq_ka.rhs().subs(vdict)\n",
    "vdict[h_c] = eq_hc.rhs().subs(vdict)\n",
    "\n",
    "\n",
    "vdict[T_l] = 305.65\n",
    "vdict[H_l] = eq_Hl.rhs().subs(vdict)\n",
    "vdict[R_ll] = eq_Rll.rhs().subs(vdict)\n",
    "\n",
    "eq_El_enbal.subs(vdict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saving definitions to separate file\n",
    "In the below, we save the definitions and variables to separate files in the /temp directory, one with the extension .sage, from which we can selectively load functions using\n",
    "`%load fun_name filenam.sage`\n",
    "and one with the extension .sobj, to be loaded elsewhere using \n",
    "`load_session()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "jupyter nbconvert  --to=python 'leaf_enbalance_eqs.ipynb'\n",
      "Exporting specified worksheet to .py file...\n",
      "Checking if specified ipynb file was run with sage kernel...\n",
      "Renaming .py file to .sage if notebook kernel was sage (to avoid exponent error)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "nbconvert returned 0\n"
     ]
    }
   ],
   "source": [
    "fun_export_ipynb('leaf_enbalance_eqs', 'temp/')\n",
    "save_session('temp/leaf_enbalance_eqs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of symbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"notruncate\">\n",
       "<table  class=\"table_form\">\n",
       "<tbody>\n",
       "<tr>\n",
       "<th>Variable</th>\n",
       "<th>Description (value)</th>\n",
       "<th>Units</th>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">a_{s}</script></td>\n",
       "<td>Fraction of one-sided leaf area covered by stomata (1 if stomata are on one side only, 2 if they are on both sides)</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{a_{sh}}</script></td>\n",
       "<td>Fraction of projected area exchanging sensible heat with the air (2)</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">\\alpha_{a}</script></td>\n",
       "<td>Thermal diffusivity of dry air</td>\n",
       "<td> m<script type=\"math/tex\">^{2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{c_{pa}}</script></td>\n",
       "<td>Specific heat of dry air (1010) </td>\n",
       "<td>J  K<script type=\"math/tex\">^{-1}</script>  kg<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{C_{wa}}</script></td>\n",
       "<td>Concentration of water in the free air </td>\n",
       "<td>mol  m<script type=\"math/tex\">^{-3}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{C_{wl}}</script></td>\n",
       "<td>Concentration of water in the leaf air space </td>\n",
       "<td>mol  m<script type=\"math/tex\">^{-3}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{D_{va}}</script></td>\n",
       "<td>Binary diffusion coefficient of water vapour in air</td>\n",
       "<td> m<script type=\"math/tex\">^{2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">E_{l}</script></td>\n",
       "<td>Latent heat flux from leaf</td>\n",
       "<td>J  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{E_{l,mol}}</script></td>\n",
       "<td>Transpiration rate in molar units</td>\n",
       "<td>mol  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">\\epsilon_{l}</script></td>\n",
       "<td>Longwave emmissivity of the leaf surface (1.0)</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">g</script></td>\n",
       "<td>Gravitational acceleration (9.81)</td>\n",
       "<td>m  s<script type=\"math/tex\">^{-2}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{g_{bw}}</script></td>\n",
       "<td>Boundary layer conductance to water vapour </td>\n",
       "<td>m  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{g_{bw,mol}}</script></td>\n",
       "<td>Boundary layer conductance to water vapour </td>\n",
       "<td>mol  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{g_{sw}}</script></td>\n",
       "<td>Stomatal conductance to water vapour</td>\n",
       "<td>m  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{g_{sw,mol}}</script></td>\n",
       "<td>Stomatal conductance to water vapour</td>\n",
       "<td>mol  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{g_{tw}}</script></td>\n",
       "<td>Total leaf conductance to water vapour</td>\n",
       "<td>m  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{g_{tw,mol}}</script></td>\n",
       "<td>Total leaf layer conductance to water vapour</td>\n",
       "<td>mol  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">h_{c}</script></td>\n",
       "<td>Average 1-sided convective transfer coefficient</td>\n",
       "<td>J  K<script type=\"math/tex\">^{-1}</script>  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">H_{l}</script></td>\n",
       "<td>Sensible heat flux from leaf</td>\n",
       "<td>J  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">k_{a}</script></td>\n",
       "<td>Thermal conductivity of dry air</td>\n",
       "<td>J  K<script type=\"math/tex\">^{-1}</script>  m<script type=\"math/tex\">^{-1}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">L_{l}</script></td>\n",
       "<td>Characteristic length scale for convection (size of leaf)</td>\n",
       "<td>m </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">\\lambda_{E}</script></td>\n",
       "<td>Latent heat of evaporation (2.45e6)</td>\n",
       "<td>J  kg<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">M_{N_{2}}</script></td>\n",
       "<td>Molar mass of nitrogen (0.028)</td>\n",
       "<td>kg  mol<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">M_{O_{2}}</script></td>\n",
       "<td>Molar mass of oxygen (0.032)</td>\n",
       "<td>kg  mol<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">M_{w}</script></td>\n",
       "<td>Molar mass of water (0.018)</td>\n",
       "<td>kg  mol<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{N_{Gr_L}}</script></td>\n",
       "<td>Grashof number</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{N_{Le}}</script></td>\n",
       "<td>Lewis number</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{N_{Nu_L}}</script></td>\n",
       "<td>Nusselt number</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{N_{Re_c}}</script></td>\n",
       "<td>Critical Reynolds number for the onset of turbulence</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{N_{Re_L}}</script></td>\n",
       "<td>Reynolds number</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{N_{Sh_L}}</script></td>\n",
       "<td>Sherwood number</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">\\nu_{a}</script></td>\n",
       "<td>Kinematic viscosity of dry air</td>\n",
       "<td> m<script type=\"math/tex\">^{2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">P_{a}</script></td>\n",
       "<td>Air pressure</td>\n",
       "<td>Pa </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{P_{N2}}</script></td>\n",
       "<td>Partial pressure of nitrogen in the atmosphere</td>\n",
       "<td>Pa </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{P_{O2}}</script></td>\n",
       "<td>Partial pressure of oxygen in the atmosphere</td>\n",
       "<td>Pa </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{P_{wa}}</script></td>\n",
       "<td>Vapour pressure in the atmosphere</td>\n",
       "<td>Pa </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{P_{was}}</script></td>\n",
       "<td>Saturation vapour pressure at air temperature</td>\n",
       "<td>Pa </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{P_{wl}}</script></td>\n",
       "<td>Vapour pressure inside the leaf</td>\n",
       "<td>Pa </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{N_{Pr}}</script></td>\n",
       "<td>Prandtl number (0.71)</td>\n",
       "<td>1 </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{r_{bw}}</script></td>\n",
       "<td>Boundary layer resistance to water vapour, inverse of <script type=\"math/tex\">g_{bw}</script></td>\n",
       "<td>s  m<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{R_{ll}}</script></td>\n",
       "<td>Longwave radiation away from leaf</td>\n",
       "<td>J  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{R_{mol}}</script></td>\n",
       "<td>Molar gas constant (8.314472)</td>\n",
       "<td>J  K<script type=\"math/tex\">^{-1}</script>  mol<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">R_{s}</script></td>\n",
       "<td>Solar shortwave flux</td>\n",
       "<td>J  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{r_{sw}}</script></td>\n",
       "<td>Stomatal resistance to water vapour, inverse of <script type=\"math/tex\">g_{sw}</script></td>\n",
       "<td>s  m<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">{r_{tw}}</script></td>\n",
       "<td>Total leaf resistance to water vapour, <script type=\"math/tex\">r_{bv} + r_{sv}</script></td>\n",
       "<td>s  m<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">\\rho_{a}</script></td>\n",
       "<td>Density of dry air</td>\n",
       "<td>kg  m<script type=\"math/tex\">^{-3}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">\\rho_{\\mathit{al}}</script></td>\n",
       "<td>Density of air at the leaf surface</td>\n",
       "<td>kg  m<script type=\"math/tex\">^{-3}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">{\\sigma}</script></td>\n",
       "<td>Stefan-Boltzmann constant (5.67e-8)</td>\n",
       "<td>J  K<script type=\"math/tex\">^{-4}</script>  m<script type=\"math/tex\">^{-2}</script>  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">T_{a}</script></td>\n",
       "<td>Air temperature</td>\n",
       "<td>K </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">T_{l}</script></td>\n",
       "<td>Leaf temperature</td>\n",
       "<td>K </td>\n",
       "</tr>\n",
       "<tr class =\"row-a\">\n",
       "<td><script type=\"math/tex\">T_{w}</script></td>\n",
       "<td>Radiative temperature of objects surrounding the leaf</td>\n",
       "<td>K </td>\n",
       "</tr>\n",
       "<tr class =\"row-b\">\n",
       "<td><script type=\"math/tex\">v_{w}</script></td>\n",
       "<td>Wind velocity</td>\n",
       "<td>m  s<script type=\"math/tex\">^{-1}</script> </td>\n",
       "</tr>\n",
       "</tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Variable    Description (value)                                                                                                   Units\n",
       "+-----------+---------------------------------------------------------------------------------------------------------------------+----------------------------------+\n",
       "  a_s         Fraction of one-sided leaf area covered by stomata (1 if stomata are on one side only, 2 if they are on both sides)   1\n",
       "  a_sh        Fraction of projected area exchanging sensible heat with the air (2)                                                  1\n",
       "  alpha_a     Thermal diffusivity of dry air                                                                                         m$^{2}$  s$^{-1}$\n",
       "  c_pa        Specific heat of dry air (1010)                                                                                       J  K$^{-1}$  kg$^{-1}$\n",
       "  C_wa        Concentration of water in the free air                                                                                mol  m$^{-3}$\n",
       "  C_wl        Concentration of water in the leaf air space                                                                          mol  m$^{-3}$\n",
       "  D_va        Binary diffusion coefficient of water vapour in air                                                                    m$^{2}$  s$^{-1}$\n",
       "  E_l         Latent heat flux from leaf                                                                                            J  m$^{-2}$  s$^{-1}$\n",
       "  E_lmol      Transpiration rate in molar units                                                                                     mol  m$^{-2}$  s$^{-1}$\n",
       "  epsilon_l   Longwave emmissivity of the leaf surface (1.0)                                                                        1\n",
       "  g           Gravitational acceleration (9.81)                                                                                     m  s$^{-2}$\n",
       "  g_bw        Boundary layer conductance to water vapour                                                                            m  s$^{-1}$\n",
       "  g_bwmol     Boundary layer conductance to water vapour                                                                            mol  m$^{-2}$  s$^{-1}$\n",
       "  g_sw        Stomatal conductance to water vapour                                                                                  m  s$^{-1}$\n",
       "  g_swmol     Stomatal conductance to water vapour                                                                                  mol  m$^{-2}$  s$^{-1}$\n",
       "  g_tw        Total leaf conductance to water vapour                                                                                m  s$^{-1}$\n",
       "  g_twmol     Total leaf layer conductance to water vapour                                                                          mol  m$^{-2}$  s$^{-1}$\n",
       "  h_c         Average 1-sided convective transfer coefficient                                                                       J  K$^{-1}$  m$^{-2}$  s$^{-1}$\n",
       "  H_l         Sensible heat flux from leaf                                                                                          J  m$^{-2}$  s$^{-1}$\n",
       "  k_a         Thermal conductivity of dry air                                                                                       J  K$^{-1}$  m$^{-1}$  s$^{-1}$\n",
       "  L_l         Characteristic length scale for convection (size of leaf)                                                             m\n",
       "  lambda_E    Latent heat of evaporation (2.45e6)                                                                                   J  kg$^{-1}$\n",
       "  M_N2        Molar mass of nitrogen (0.028)                                                                                        kg  mol$^{-1}$\n",
       "  M_O2        Molar mass of oxygen (0.032)                                                                                          kg  mol$^{-1}$\n",
       "  M_w         Molar mass of water (0.018)                                                                                           kg  mol$^{-1}$\n",
       "  Gr          Grashof number                                                                                                        1\n",
       "  Le          Lewis number                                                                                                          1\n",
       "  Nu          Nusselt number                                                                                                        1\n",
       "  Re_c        Critical Reynolds number for the onset of turbulence                                                                  1\n",
       "  Re          Reynolds number                                                                                                       1\n",
       "  Sh          Sherwood number                                                                                                       1\n",
       "  nu_a        Kinematic viscosity of dry air                                                                                         m$^{2}$  s$^{-1}$\n",
       "  P_a         Air pressure                                                                                                          Pa\n",
       "  P_N2        Partial pressure of nitrogen in the atmosphere                                                                        Pa\n",
       "  P_O2        Partial pressure of oxygen in the atmosphere                                                                          Pa\n",
       "  P_wa        Vapour pressure in the atmosphere                                                                                     Pa\n",
       "  P_was       Saturation vapour pressure at air temperature                                                                         Pa\n",
       "  P_wl        Vapour pressure inside the leaf                                                                                       Pa\n",
       "  Pr          Prandtl number (0.71)                                                                                                 1\n",
       "  r_bw        Boundary layer resistance to water vapour, inverse of $g_{bw}$                                                        s  m$^{-1}$\n",
       "  R_ll        Longwave radiation away from leaf                                                                                     J  m$^{-2}$  s$^{-1}$\n",
       "  R_mol       Molar gas constant (8.314472)                                                                                         J  K$^{-1}$  mol$^{-1}$\n",
       "  R_s         Solar shortwave flux                                                                                                  J  m$^{-2}$  s$^{-1}$\n",
       "  r_sw        Stomatal resistance to water vapour, inverse of $g_{sw}$                                                              s  m$^{-1}$\n",
       "  r_tw        Total leaf resistance to water vapour, $r_{bv} + r_{sv}$                                                              s  m$^{-1}$\n",
       "  rho_a       Density of dry air                                                                                                    kg  m$^{-3}$\n",
       "  rho_al      Density of air at the leaf surface                                                                                    kg  m$^{-3}$\n",
       "  sigm        Stefan-Boltzmann constant (5.67e-8)                                                                                   J  K$^{-4}$  m$^{-2}$  s$^{-1}$\n",
       "  T_a         Air temperature                                                                                                       K\n",
       "  T_l         Leaf temperature                                                                                                      K\n",
       "  T_w         Radiative temperature of objects surrounding the leaf                                                                 K\n",
       "  v_w         Wind velocity                                                                                                         m  s$^{-1}$"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Creating dictionary to substitute names of units with shorter forms\n",
    "var('m s J Pa K kg mol')\n",
    "subsdict = {meter: m, second: s, joule: J, pascal: Pa, kelvin: K, kilogram: kg, mole: mol}\n",
    "var('N_Re_L N_Re_c N_Le N_Nu_L N_Gr_L N_Sh_L')\n",
    "dict_varnew = {Re: N_Re_L, Re_c: N_Re_c, Le: N_Le, Nu: N_Nu_L, Gr: N_Gr_L, Sh: N_Sh_L}\n",
    "dict_varold = {v: k for k, v in dict_varnew.iteritems()}\n",
    "variables = sorted([str(variable.subs(dict_varnew)) for variable in udict.keys()],key=str.lower)\n",
    "tableheader = [('Variable', 'Description (value)', 'Units')]\n",
    "tabledata = [('Variable', 'Description (value)', 'Units')]\n",
    "for variable1 in variables:\n",
    "    variable2 = eval(variable1).subs(dict_varold)\n",
    "    variable = str(variable2)\n",
    "    tabledata.append((eval(variable),docdict[eval(variable)],fun_units_formatted(variable)))\n",
    "\n",
    "table(tabledata, header_row=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 7.3",
   "language": "",
   "name": "sagemath"
  },
  "language": "python",
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}