"investigating the stress factor missfunctionning in simpler_model.ipynb"

notebooks/theory/pyFile_storage/theory_equation.py

-# -*- coding: utf-8 -*-
-#
-# This file is for use with essm.
-# Copyright (C) 2021 ETH Zurich, Swiss Data Science Center.
-#
-# essm is free software; you can redistribute it
-# and/or modify it under the terms of the GNU General Public License as
-# published by the Free Software Foundation; either version 2 of the
-# License, or (at your option) any later version.
-#
-# essm is distributed in the hope that it will be
-# useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-# General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with essm; if not, write to the
-# Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
-# MA 02111-1307, USA.
-"""Equations defined in notebooks/theory/theory.ipynb and dependencies."""
-
-from __future__ import division
-
-from __main__ import (
-    T0, VPD, Delta, E, G, L, Mw, Pvs, R, Rn, S, S_mvg, T, alpha, alpha_PT, c1,
-    c_p, e_a, e_s, g_a, g_s, gamma, h, m, n, r_a, r_s, rho_a, theta, theta_1,
-    theta_2, theta_3, theta_4, theta_res, theta_sat
-)
-from essm import Eq
-from essm.equations import Equation
-from sympy import Abs, Eq, Piecewise, exp, log
-
-
-class eq_m_n(Equation):
-    """Relation between n and m parameters in MVG model"""
-
-    expr = Eq(m, 1 - 1 / n)
-
-
-class eq_MVG_neg_case(Equation):
-    """Mualem-Van Genuchtem model if h < 0"""
-
-    expr = Eq(
-        theta,
-        theta_res + ((-theta_res + theta_sat) / (Abs(alpha * h) ** n + 1)) ** m
-    )
-
-
-class eq_MVG(Equation):
-    """Complet Van Genuchtem model"""
-
-    expr = Eq(
-        theta,
-        Piecewise((theta_sat, h > 0), (
-            theta_res + ((-theta_res + theta_sat) /
-                         (Abs(alpha * h) ** n + 1)) ** m, True
-        ))
-    )
-
-
-class eq_sat_degree(Equation):
-    """Relative Saturation degree equation"""
-
-    expr = Eq(S_mvg, (theta - theta_res) / (-theta_res + theta_sat))
-
-
-class eq_MVG_h(Equation):
-    """Mualem Van Genuchten model, h as function of theta"""
-
-    expr = Eq(h, (-1 + S_mvg ** (-1 / m)) ** (1 / n) / alpha)
-
-
-class eq_h_FC(Equation):
-    """Compute the field capacity out of the soil properties"""
-
-    expr = Eq(h, ((n - 1) / n) ** ((1 - 2 * n) / n) / alpha)
-
-
-class eq_theta_4_3(Equation):
-    """Equation valide between theta_4 and theta_3"""
-
-    expr = Eq(
-        theta, theta / (theta_3 - theta_4) + theta_4 / (-theta_3 + theta_4)
-    )
-
-
-class eq_theta_2_1(Equation):
-    """Equation valide between theta_2 and theta_1"""
-
-    expr = Eq(
-        theta, theta / (-theta_1 + theta_2) + theta_1 / (theta_1 - theta_2)
-    )
-
-
-class eq_water_stress_simple(Equation):
-    """Simple water stress factor function (h_3 independant of T_0)"""
-
-    expr = Eq(
-        S,
-        Piecewise((0, theta < theta_4), (
-            theta / (theta_3 - theta_4) + theta_4 /
-            (-theta_3 + theta_4), theta < theta_3
-        ), (1, theta < theta_2), (
-            theta / (-theta_1 + theta_2) + theta_1 /
-            (theta_1 - theta_2), theta < theta_1
-        ), (0, True))
-    )
-
-
-class eq_Pvs_T(Equation):
-    """Saturation vapour pressure, Slide 8"""
-
-    expr = Eq(Pvs, c1 * exp(L * Mw * (1 / T0 - 1 / T) / R))
-
-
-class eq_Delta(Equation):
-    """Slope of d.Pvas(T)/d.T """
-
-    expr = Eq(
-        Delta,
-        L * Mw * c1 * exp(L * Mw * (1 / T0 - 1 / T) / R) / (R * T ** 2)
-    )
-
-
-class eq_PT(Equation):
-    """Priestley-Taylor equation for evaporation flux"""
-
-    expr = Eq(E, Delta * alpha_PT * (-G + Rn) / (L * (Delta + gamma)))
-
-
-class eq_PM(Equation):
-    """Penman-Monteith equation"""
-
-    expr = Eq(
-        E, (Delta * (-G + Rn) + c_p * rho_a * (-e_a + e_s) / r_a) /
-        (L * (Delta + gamma * (1 + r_s / r_a)))
-    )
-
-
-class eq_PM_VPD(Equation):
-    """Penman-Monteith equation using VPD"""
-
-    expr = Eq(
-        E, (Delta * (-G + Rn) + VPD * c_p * rho_a / r_a) /
-        (L * (Delta + gamma * (1 + r_s / r_a)))
-    )
-
-
-class eq_PM_g(Equation):
-    """Penman-Monteith equation using stomatal conductance"""
-
-    expr = Eq(
-        E, (Delta * (-G + Rn) + VPD * c_p * g_a * rho_a) /
-        (L * (Delta + gamma * (g_a / g_s + 1)))
-    )
-
-
-class eq_PM_inv(Equation):
-    """Inverse PM equation for inverse modelling of surface resistance"""
-
-    expr = Eq(
-        r_s, (
-            -Delta * E * L * r_a - Delta * G * r_a + Delta * Rn * r_a -
-            E * L * gamma * r_a + VPD * c_p * rho_a
-        ) / (E * L * gamma)
-    )
-
-
-__all__ = (
-    'eq_m_n',
-    'eq_MVG_neg_case',
-    'eq_MVG',
-    'eq_sat_degree',
-    'eq_MVG_h',
-    'eq_h_FC',
-    'eq_theta_4_3',
-    'eq_theta_2_1',
-    'eq_water_stress_simple',
-    'eq_Pvs_T',
-    'eq_Delta',
-    'eq_PT',
-    'eq_PM',
-    'eq_PM_VPD',
-    'eq_PM_g',
-    'eq_PM_inv',
-)