Tutorial 1: Basic model run

from fastcore import *
from nbdev.showdoc import *

from bonan_hydraulics.utils import satvap
from bonan_hydraulics.root_parms import root_params
from bonan_hydraulics.soil_params import soil_params
from bonan_hydraulics.leaf_fluxes import leaf_fluxes
from bonan_hydraulics.soil_resistance import soil_resistance
from bonan_hydraulics.leaf_phys_params import leaf_phys_params
from bonan_hydraulics.parameter_classes import (
    PhysCon,
    Params,
    Ground,
    Soil,
    RootVar,
    Leaf,
    Flux,
    Atmos,
)

def run_bonan_hydraulics_model():
    """Run the plant hydraulics model (Supplemental Program 13.1)."""

    # Waveband indices ----------------------------------------------------------
    params = Params(vis=0, nir=1)

    # Physical constants --------------------------------------------------------
    physcon = PhysCon()

    # Atmospheric forcing -------------------------------------------------------
    leaftype = "sun"
    atmos = Atmos()
    atmos.co2air = 380.0
    atmos.o2air = 0.209 * 1000.0

    atmos.tair = physcon.tfrz + 30.0
    atmos.relhum = 60.0

    atmos.wind = 1.0
    atmos.patm = 101325.0

    # Vapor pressure and specific humidity
    esat, desat = satvap(atmos.tair - physcon.tfrz)
    atmos.eair = esat * (atmos.relhum / 100.0)
    atmos.qair = (
        physcon.mmh2o
        / physcon.mmdry
        * atmos.eair
        / (atmos.patm - (1.0 - physcon.mmh2o / physcon.mmdry) * atmos.eair)
    )

    # Ground variables ----------------------------------------------------------
    ground = Ground()
    ground.albsoi[params.vis] = 0.1
    ground.albsoi[params.nir] = 0.2
    tg = atmos.tair
    ground.irgrd = physcon.sigma * tg**4

    # Set soil ------------------------------------------------------------------
    soil = Soil()

    # Loam
    soil.texture = 5
    soil = soil_params(soil)

    # Set soil moisture ---------------------------------------------------------
    soil.h2osoi_vol = []
    soil.psi = []
    for j in range(soil.nlevsoi):
        h2osoi = 0.5 * soil.watsat[j]
        soil.h2osoi_vol.append(h2osoi)
        soil.psi.append(soil.psisat[j] * (h2osoi / soil.watsat[j]) ** (-soil.bsw[j]))

    # Set root parameters -------------------------------------------------------
    rootvar = root_params()

    # Set site variables --------------------------------------------------------
    flux = Flux()
    flux.height = 15.0  # Leaf height (m)
    rootvar.biomass = 500.0  # Fine root biomass (g biomass / m2)
    flux.lai = 5.0  # Canopy leaf area index (m2/m2)

    # Set leaf physiology variables ---------------------------------------------
    leaf = Leaf()

    # C3 plant
    leaf.c3psn = 1

    # Co-limitation on
    leaf.colim = 1

    leaf = leaf_phys_params(params, physcon, leaf)

    # Radiation absorbed by leaf ------------------------------------------------
    flux.swinc[params.vis] = atmos.swsky[params.vis] * (1.0 + ground.albsoi[params.vis])
    flux.swinc[params.nir] = atmos.swsky[params.nir] * (1.0 + ground.albsoi[params.nir])

    flux.swflx[params.vis] = flux.swinc[params.vis] * (
        1.0 - leaf.rho[params.vis] - leaf.tau[params.vis]
    )
    flux.swflx[params.nir] = flux.swinc[params.nir] * (
        1.0 - leaf.rho[params.nir] - leaf.tau[params.nir]
    )
    flux.apar = flux.swflx[params.vis] * 4.6

    flux.qa = (
        flux.swflx[params.vis]
        + flux.swflx[params.nir]
        + leaf.emiss * (atmos.irsky + ground.irgrd)
    )

    # Hydraulic conductances ----------------------------------------------------
    flux.rplant = 1.0 / leaf.gplant
    flux = soil_resistance(physcon, leaf, rootvar, soil, flux)
    flux.lsc = 1.0 / (flux.rsoil + flux.rplant)

    print(f"\n--- Hydraulic conductances ---")
    print(f"rplant = {flux.rplant:15.5f}")
    print(f"rsoil  = {flux.rsoil:15.5f}")
    print(f"lsc    = {flux.lsc:15.5f}")

    # Molar density (mol/m3)
    atmos.rhomol = atmos.patm / (physcon.rgas * atmos.tair)

    # Air density (kg/m3)
    atmos.rhoair = (
        atmos.rhomol
        * physcon.mmdry
        * (1.0 - (1.0 - physcon.mmh2o / physcon.mmdry) * atmos.eair / atmos.patm)
    )

    # Molecular mass of air (kg/mol)
    atmos.mmair = atmos.rhoair / atmos.rhomol

    # Specific heat of air at constant pressure (J/mol/K)
    atmos.cpair = (
        physcon.cpd
        * (1.0 + (physcon.cpw / physcon.cpd - 1.0) * atmos.qair)
        * atmos.mmair
    )

    # Atmospheric longwave radiation
    atmos.irsky = 400.0

    # Solar radiation
    if leaftype == "sun":
        fsds = 800.0
    else:
        fsds = 300.0
    atmos.swsky[params.vis] = 0.5 * fsds
    atmos.swsky[params.nir] = 0.5 * fsds

    # Initial conditions --------------------------------------------------------
    flux.tleaf = atmos.tair
    flux.psi_leaf = -0.1

    # 30 minutes in seconds
    flux.dt = 30.0 * 60.0

    # Calculate fluxes ----------------------------------------------------------
    flux = leaf_fluxes(physcon, atmos, leaf, flux)

    # Print output --------------------------------------------------------------
    print(f"\n----------- Results -----------")
    print(f"dleaf    = {leaf.dleaf * 100:15.5f} cm")
    print(f"apar     = {flux.apar:15.5f} umol/m2/s")
    print(f"tleaf    = {flux.tleaf - physcon.tfrz:15.5f} degC")
    print(f"qa       = {flux.qa:15.5f} W/m2")
    print(f"lhflx    = {flux.lhflx:15.5f} W/m2")
    print(f"etflx    = {flux.etflx * 1000:15.5f} mmol H2O/m2/s")
    print(f"an       = {flux.an:15.5f} umol CO2/m2/s")
    print(f"gbh      = {flux.gbh:15.5f} mol/m2/s")
    print(f"gs       = {flux.gs:15.5f} mol H2O/m2/s")
    print(f"psi_leaf = {flux.psi_leaf:15.5f} MPa")


if __name__ == "__main__":
    run_bonan_hydraulics_model()

--- Hydraulic conductances ---
rplant =         0.25000
rsoil  =         0.25012
lsc    =         1.99952

----------- Results -----------
dleaf    =         5.00000 cm
apar     =         0.00000 umol/m2/s
tleaf    =        22.74317 degC
qa       =       469.28797 W/m2
lhflx    =         0.18930 W/m2
etflx    =         0.00432 mmol H2O/m2/s
an       =        -0.83307 umol CO2/m2/s
gbh      =         0.90262 mol/m2/s
gs       =         0.00200 mol H2O/m2/s
psi_leaf =        -0.28756 MPa