R/RcppExports.R

In medfate: Mediterranean Forest Simulation

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Physical and biophysical utility functions
#' 
#' Utility functions for the calculation of biophysical variables. 
#'
#' @param t Time of the day (in seconds).
#' @param daylength Day length (in seconds).
#' 
#' @details 
#' Functions \code{biophysics_leafTemperature} and \code{biophysics_leafTemperature2} calculate leaf temperature according to energy balance equation given in Campbell and Norman (1988). 
#' 
#' Function \code{biophysics_radiationDiurnalPattern} follows the equations given in Liu and Jordan (1960). 
#' 
#' Function \code{biophysics_temperatureDiurnalPattern} determines diurnal temperature pattern assuming a sinusoidal pattern with T = Tmin at sunrise and T = (Tmin+Tmax)/2 at sunset and a linear change in temperature between sunset and Tmin of the day after (McMurtrie et al. 1990). 
#' 
#' Function \code{biophysics_waterDynamicViscosity} calculates water dynamic viscosity following the Vogel (1921) equation.
#' 
#' @return
#' Values returned for each function are:
#' \itemize{
#'   \item{\code{biophysics_leafTemperature} and \code{biophysics_leafTemperature2}: leaf temperature (in ºC)}
#'   \item{\code{biophysics_leafVapourPressure}: leaf vapour pressure (in kPa)} 
#'   \item{\code{biophysics_radiationDiurnalPattern}: the proportion of daily radiation corresponding to the input time in seconds after sunrise.} 
#'   \item{\code{biophysics_temperatureDiurnalPattern}: diurnal pattern of temperature.}
#'   \item{\code{biophysics_waterDynamicViscosity}: Water dynamic viscosity relative to 20ºC.} 
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @references
#' Campbell, G. S., and J. M. Norman. 1998. An introduction to environmental biophysics: 2nd edition. (eqns. 14.1 & 14.3)
#'   
#' B. Y. H. Liu and R. C. Jordan, “The interrelationship and characteristic distribution of direct, diffuse and total solar radiation,” Solar Energy, vol. 4, no. 3, pp. 1–19, 1960. 
#' 
#' McMurtrie, R. E., D. A. Rook, and F. M. Kelliher. 1990. Modelling the yield of Pinus radiata on a site limited by water and nitrogen. Forest Ecology and Management 30:381–413.
#' 
#' H. Vogel, "Das Temperaturabhangigkeitsgesetz der Viskositat von Flussigkeiten", Physikalische Zeitschrift, vol. 22, pp. 645–646, 1921.
#' 
#' @seealso \code{\link{spwb}}
#' 
#' @name biophysics
biophysics_radiationDiurnalPattern <- function(t, daylength) {
    .Call(`_medfate_radiationDiurnalPattern`, t, daylength)
}

#' @rdname biophysics
#' @param tmin,tmax Minimum and maximum daily temperature (ºC).
#' @param tminPrev,tmaxPrev,tminNext Maximum and minimum daily temperatures of the previous and following day (ºC).
biophysics_temperatureDiurnalPattern <- function(t, tmin, tmax, tminPrev, tmaxPrev, tminNext, daylength) {
    .Call(`_medfate_temperatureDiurnalPattern`, t, tmin, tmax, tminPrev, tmaxPrev, tminNext, daylength)
}

#' @rdname biophysics
#' 
#' @param u Wind speed above the leaf boundary layer (in m/s).
#' @param airTemperature Air temperature (in ºC).
#' @param absRad Absorbed long- and short-wave radiation (in W·m-2).
#' @param E Transpiration flow (in mmol H20·m-2·s-1) per one sided leaf area basis.
#' @param leafWidth Leaf width (in cm).
#' 
biophysics_leafTemperature <- function(absRad, airTemperature, u, E, leafWidth = 1.0) {
    .Call(`_medfate_leafTemperature`, absRad, airTemperature, u, E, leafWidth)
}

#' @rdname biophysics
#' @param SWRabs Absorbed short-wave radiation (in W·m-2).
#' @param LWRnet Net long-wave radiation balance (in W·m-2).
biophysics_leafTemperature2 <- function(SWRabs, LWRnet, airTemperature, u, E, leafWidth = 1.0) {
    .Call(`_medfate_leafTemperature2`, SWRabs, LWRnet, airTemperature, u, E, leafWidth)
}

#' @rdname biophysics
#' @param leafTemp Leaf temperature (ºC).
#' @param leafPsi Leaf water potential (MPa).
biophysics_leafVapourPressure <- function(leafTemp, leafPsi) {
    .Call(`_medfate_leafVapourPressure`, leafTemp, leafPsi)
}

#' @rdname biophysics
#' @param I Irradiance (in W*m-2).
#' @param lambda Wavelength (in nm).
biophysics_irradianceToPhotonFlux <- function(I, lambda = 546.6507) {
    .Call(`_medfate_irradianceToPhotonFlux`, I, lambda)
}

#' @rdname biophysics
#' @param temp Temperature (ºC).
biophysics_waterDynamicViscosity <- function(temp) {
    .Call(`_medfate_waterDynamicViscosity`, temp)
}

#' Carbon-related functions
#' 
#' Low-level functions used in the calculation of carbon balance.
#' 
#' @param LAI Leaf area index.
#' @param N Density (ind·ha-1).
#' @param SLA Specific leaf area (mm2/mg = m2/kg).
#' @param leafDensity  Density of leaf tissue (dry weight over volume).
#' @param SA Sapwood area (cm2).
#' @param H Plant height (cm).
#' @param L Coarse root length (mm) for each soil layer.
#' @param V Proportion of fine roots in each soil layer.
#' @param woodDensity Wood density (dry weight over volume).
#' @param conduit2sapwood Proportion of sapwood corresponding to conducive elements (vessels or tracheids) as opposed to parenchymatic tissue.
#' @param osmoticWP Osmotic water potential (MPa).
#' @param temp Temperature (degrees Celsius).
#' @param nonSugarConc Concentration of inorganic solutes (mol/l).
#' @param sugarConc Concentration of soluble sugars (mol/l).
#' @param starchConc Concentration of starch (mol/l)
#' @param eqSugarConc Equilibrium concentration of soluble sugars (mol/l).
#' 
#' @return Values returned for each function are:
#' \itemize{
#'   \item{\code{carbon_leafStarchCapacity}: Capacity of storing starch in the leaf compartment (mol gluc/ind.).}
#'   \item{\code{carbon_leafStructuralBiomass}: Leaf structural biomass (g dry/ind.)}
#'   \item{\code{carbon_sapwoodStarchCapacity}: Capacity of storing starch in the sapwood compartment (mol gluc/ind.).}
#'   \item{\code{carbon_sapwoodStructuralBiomass}: Sapwood structural biomass (g dry/ind.)}
#'   \item{\code{carbon_sapwoodStructuralLivingBiomass}: Living sapwood (parenchyma) structural biomass (g dry/ind.)}
#'   \item{\code{carbon_sugarConcentration}: Sugar concentration (mol gluc/l)}
#'   \item{\code{carbon_osmoticWaterPotential}: Osmotic component of water potential (MPa)}
#'   \item{\code{carbon_relativeSapViscosity}: Relative viscosity of sapwood with respect to pure water (according to Forst et al. (2002)).}
#'   \item{\code{carbon_sugarStarchDynamicsLeaf}: Rate of conversion from sugar to starch in leaf (mol gluc/l/s).}
#'   \item{\code{carbon_sugarStarchDynamicsStem}: Rate of conversion from sugar to starch in leaf (mol gluc/l/s).}
#'   \item{\code{carbon_carbonCompartments}: A data frame with the size of compartments for each plant cohort, in the specified units.}
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @references
#' Forst P, Wermer F, Delgado A (2002). On the pressure dependence of the viscosity of aqueous sugar solutions. Rheol Acta 41: 369–374 DOI 10.1007/s00397-002-0238-y
#' 
#' @seealso \code{\link{growth}}
#' 
#' @name carbon
carbon_sugarStarchDynamicsLeaf <- function(sugarConc, starchConc, eqSugarConc) {
    .Call(`_medfate_sugarStarchDynamicsLeaf`, sugarConc, starchConc, eqSugarConc)
}

#' @rdname carbon
carbon_sugarStarchDynamicsStem <- function(sugarConc, starchConc, eqSugarConc) {
    .Call(`_medfate_sugarStarchDynamicsStem`, sugarConc, starchConc, eqSugarConc)
}

#' @rdname carbon
carbon_osmoticWaterPotential <- function(sugarConc, temp, nonSugarConc) {
    .Call(`_medfate_osmoticWaterPotential`, sugarConc, temp, nonSugarConc)
}

#' @rdname carbon
carbon_sugarConcentration <- function(osmoticWP, temp, nonSugarConc) {
    .Call(`_medfate_sugarConcentration`, osmoticWP, temp, nonSugarConc)
}

#' @rdname carbon
carbon_relativeSapViscosity <- function(sugarConc, temp) {
    .Call(`_medfate_relativeSapViscosity`, sugarConc, temp)
}

#' @rdname carbon
carbon_leafStructuralBiomass <- function(LAI, N, SLA) {
    .Call(`_medfate_leafStructuralBiomass`, LAI, N, SLA)
}

#' @rdname carbon
carbon_leafStarchCapacity <- function(LAI, N, SLA, leafDensity) {
    .Call(`_medfate_leafStarchCapacity`, LAI, N, SLA, leafDensity)
}

#' @rdname carbon
carbon_sapwoodStructuralBiomass <- function(SA, H, L, V, woodDensity) {
    .Call(`_medfate_sapwoodStructuralBiomass`, SA, H, L, V, woodDensity)
}

#' @rdname carbon
carbon_sapwoodStructuralLivingBiomass <- function(SA, H, L, V, woodDensity, conduit2sapwood) {
    .Call(`_medfate_sapwoodStructuralLivingBiomass`, SA, H, L, V, woodDensity, conduit2sapwood)
}

#' @rdname carbon
carbon_sapwoodStarchCapacity <- function(SA, H, L, V, woodDensity, conduit2sapwood) {
    .Call(`_medfate_sapwoodStarchCapacity`, SA, H, L, V, woodDensity, conduit2sapwood)
}

#' @rdname carbon
#' @param x An object of class \code{\link{growthInput}}.
#' @param biomassUnits A string for output biomass units, either "g_ind" (g per individual) or "g_m2" (g per square meter).
carbon_carbonCompartments <- function(x, biomassUnits = "g_m2") {
    .Call(`_medfate_carbonCompartments`, x, biomassUnits)
}

.criticalFirelineIntensity <- function(CBH, M) {
    .Call(`_medfate_criticalFirelineIntensity`, CBH, M)
}

#' Fire behaviour functions
#' 
#' Function \code{fire_FCCS()} implements a modification of the fire behavior models 
#' described for the Fuel Characteristics Classification System (FCCS) in Prichard et al. (2013). 
#' Function \code{fire_Rothermel()} implements Rothermel's (1972) fire behaviour 
#' model (modified from package 'Rothermel' (Giorgio Vacchiano, Davide Ascoli)).
#' 
#' @param FCCSpropsSI A data frame describing the properties of five fuel strata (canopy, shrub, herbs, dead woody and litter) returned by \code{\link{fuel_FCCS}}.
#' @param MliveSI Moisture of live fuels (in percent of dry weight) for canopy, shrub, and herb strata. Live moisture values are drawn from column \code{ActFCM} in \code{FCCSpropsSI} if available (see \code{\link{fuel_FCCS}}). Otherwise, moisture values supplied for \code{MliveSI} are used.
#' @param MdeadSI Moisture of dead fuels (in percent of dry weight) for canopy, shrub, herb, woody and litter strata.
#' @param slope Slope (in degrees).
#' @param windSpeedSI Wind speed (in m/s) at 20 ft (6 m) over vegetation (default 11 m/s = 40 km/h)
#' 
#' @details Default moisture, slope and windspeed values are benchmark conditions 
#' used to calculate fire potentials (Sandberg et al. 2007) and map vulnerability to fire.
#' 
#' @return Both functions return list with fire behavior variables. 
#' 
#' In the case of \code{fire_FCCS}, the function returns the variables in three blocks (lists \code{SurfaceFire}, \code{CrownFire} and \code{FirePotentials}), and the values are:
#' \itemize{
#'   \item{\code{SurfaceFire$`midflame_WindSpeed [m/s]`}: Midflame wind speed in the surface fire.}
#'   \item{\code{SurfaceFire$phi_wind}: Spread rate modifier due to wind.}
#'   \item{\code{SurfaceFire$phi_slope}: Spread rate modifier due to slope.}
#'   \item{\code{SurfaceFire$`I_R_surf [kJ/m2/min]`}: Intensity of the surface fire reaction.}
#'   \item{\code{SurfaceFire$`I_R_litter [kJ/m2/min]`}: Intensity of the litter fire reaction.}
#'   \item{\code{SurfaceFire$`q_surf [kJ/m2]`}: Heat sink of the surface fire.}
#'   \item{\code{SurfaceFire$`q_litter [kJ/m2]`}: Heat sink of the litter fire.}
#'   \item{\code{SurfaceFire$xi_surf}: Propagating flux ratio of the surface fire.}
#'   \item{\code{SurfaceFire$xi_litter}: Propagating flux ratio of the litter fire.}
#'   \item{\code{SurfaceFire$`ROS_surf [m/min]`}: Spread rate of the surface fire(without accounting for faster spread in the litter layer).}
#'   \item{\code{SurfaceFire$`ROS_litter [m/min]`}: Spread rate of the litter fire.}
#'   \item{\code{SurfaceFire$`ROS_windslopecap [m/min]`}: Maximum surface fire spread rate according to wind speed.}
#'   \item{\code{SurfaceFire$`ROS [m/min]`}: Final spread rate of the surface fire.}
#'   \item{\code{SurfaceFire$`I_b [kW/m]`}: Fireline intensity of the surface fire.}
#'   \item{\code{SurfaceFire$`FL [m]`}: Flame length of the surface fire.}
#'   \item{\code{CrownFire$`I_R_canopy [kJ/m2/min]`}: Intensity of the canopy fire reaction.}
#'   \item{\code{CrownFire$`I_R_crown [kJ/m2/min]`}: Intensity of the crown fire reaction (adding surface and canopy reactions).}
#'   \item{\code{CrownFire$`q_canopy [kJ/m2]`}: Heat sink of the canopy fire.}
#'   \item{\code{CrownFire$`q_crown [kJ/m2]`}: Heat sink of the crown fire (adding surface and canopy heat sinks).}
#'   \item{\code{CrownFire$xi_surf}: Propagating flux ratio of the crown fire.}
#'   \item{\code{CrownFire$`canopy_WindSpeed [m/s]`}: Wind speed in the canopy fire (canopy top wind speed).}
#'   \item{\code{CrownFire$WAF}: Wind speed adjustment factor for crown fires.}
#'   \item{\code{CrownFire$`ROS [m/min]`}: Spread rate of the crown fire.}
#'   \item{\code{CrownFire$Ic_ratio}: Crown initiation ratio.}
#'   \item{\code{CrownFire$`I_b [kW/m]`}: Fireline intensity of the crown fire.}
#'   \item{\code{CrownFire$`FL [m]`}: Flame length of the crown fire.}
#'   \item{\code{FirePotentials$RP}: Surface fire reaction potential ([0-9]).}
#'   \item{\code{FirePotentials$SP}: Surface fire spread rate potential ([0-9]).}
#'   \item{\code{FirePotentials$FP}: Surface fire flame length potential ([0-9]).}
#'   \item{\code{FirePotentials$SFP}: Surface fire potential ([0-9]).}
#'   \item{\code{FirePotentials$IC}: Crown initiation potential ([0-9]).}
#'   \item{\code{FirePotentials$TC}: Crown-to-crown transmission potential ([0-9]).}
#'   \item{\code{FirePotentials$RC}: Crown fire spread rate potential ([0-9]).}
#'   \item{\code{FirePotentials$CFC}: Crown fire potential ([0-9]).}
#' }
#' 
#' @references
#' Albini, F. A. (1976). Computer-based models of wildland fire behavior: A users' manual. Ogden, UT: US Department of Agriculture, Forest Service, Intermountain Forest and Range Experiment Station.
#' 
#' Rothermel, R. C. 1972. A mathematical model for predicting fire spread in wildland fuels. USDA Forest Service Research Paper INT USA.
#' 
#' Prichard, S. J., D. V Sandberg, R. D. Ottmar, E. Eberhardt, A. Andreu, P. Eagle, and K. Swedin. 2013. Classification System Version 3.0: Technical Documentation.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @note Default moisture, slope and windspeed values are benchmark conditions used to calculate fire potentials (Sandberg et al. 2007) and map vulnerability to fire.
#' 
#' @seealso \code{\link{fuel_FCCS}}
#' 
#' @examples
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Calculate fuel properties according to FCCS
#' fccs = fuel_FCCS(exampleforestMED, SpParamsMED)
#'   
#' #Calculate fire behavior according to FCCS
#' fire_FCCS(fccs)
#'   
#' #Load fuel model parameter data
#' data(SFM_metric)
#'       
#' #Fuel stratification (returns heights in cm)
#' fs = fuel_stratification(exampleforestMED, SpParamsMED)
#' 
#' #Correct windspeed (transform heights to m)
#' u = 11 #m/s
#' umf = u*fuel_windAdjustmentFactor(fs$surfaceLayerTopHeight/100, 
#'                                   fs$canopyBaseHeight/100, 
#'                                   fs$canopyTopHeight/100, 60)
#'       
#' #Call Rothermel function using fuel model 'A6'
#' fire_Rothermel(modeltype="D", wSI = as.numeric(SFM_metric["A6",2:6]), 
#'                sSI = as.numeric(SFM_metric["A6",7:11]), 
#'                delta = as.numeric(SFM_metric["A6",12]),
#'                mx_dead = as.numeric(SFM_metric["A6",13]),
#'                hSI = as.numeric(SFM_metric["A6",14:18]),
#'                mSI = c(10,10,10,30,60),
#'                u=umf, windDir=0, slope=0, aspect=0)
#'             
#'  
#' @name fire_behaviour
fire_FCCS <- function(FCCSpropsSI, MliveSI = as.numeric( c(90, 90, 60)), MdeadSI = as.numeric( c(6, 6, 6, 6, 6)), slope = 0.0, windSpeedSI = 11.0) {
    .Call(`_medfate_FCCSbehaviour`, FCCSpropsSI, MliveSI, MdeadSI, slope, windSpeedSI)
}

#' @rdname fire_behaviour
#' 
#' @param modeltype 'S'(tatic) or 'D'(ynamic)
#' @param wSI A vector of fuel load (t/ha) for five fuel classes.
#' @param sSI A vector of surface-to-volume ratio (m2/m3) for five fuel classes.
#' @param delta A value of fuel bed depth (cm).
#' @param mx_dead A value of dead fuel moisture of extinction (percent).
#' @param hSI A vector of heat content (kJ/kg) for five fuel classes.
#' @param mSI A vector of percent moisture on a dry weight basis (percent) for five fuel classes.
#' @param u A value of windspeed (m/s) at midflame height.
#' @param windDir Wind direction (in degrees from north). North means blowing from north to south.
#' @param aspect Aspect (in degrees from north).
#' 
fire_Rothermel <- function(modeltype, wSI, sSI, delta, mx_dead, hSI, mSI, u, windDir, slope, aspect) {
    .Call(`_medfate_rothermel`, modeltype, wSI, sSI, delta, mx_dead, hSI, mSI, u, windDir, slope, aspect)
}

#' Fire severity functions
#' 
#' Functions to estimate fire effects on foliage, buds and cambium, based on the model
#' by Michaletz & Johnson (2008)
#' 
#' @param Ib_surf Surface fireline intensity (kW/m).
#' @param t_res fire residence time (seconds).
#' @param T_air Air temperature (degrees Celsius).
#' @param T_necrosis Temperature of tissue necrosis (degrees Celsius).
#' @param rho_air Air density (kg/m3).
#' @param rho_bark Bark density (kg/m3).
#' @param fmc_bark Bark moisture content (\% dry weight).
#' @param z height (m).
#' @param SLA Specific leaf area (m2/kg).
#' @param h Heat transfer coefficient
#' @param c Specific heat capacity
#' @param thermal_factor Tissue thermal factor.
#' @param bark_diffusivity Bark thermal diffusivity (m2/s).
#' 
#' @return 
#' \itemize{
#'   \item{Function \code{fire_plumeTemperature} returns the plume temperature at a given height.}
#'   \item{Function \code{fire_barkThermalDiffusivity} returns the bark thermal diffusivity given a bark moisture value.}
#'   \item{Function \code{fire_radialBoleNecrosis} returns the depth of radial bole necrosis in cm.}
#'   \item{Function \code{fire_leafThermalFactor} returns the thermal factor of leaves as a function of specific leaf area.}
#'   \item{Function \code{fire_necrosisCriticalTemperature} returns the (plume) temperature yielding necrosis for a given residence time and tissue thermal factor.}
#'   \item{Function \code{fire_necrosisHeight} returns the height (in m) of necrosis for tissues with given thermal factor.}
#' }
#' 
#' @references
#' 
#'   Michaletz, S.T., and Johnson, E.A. 2006. A heat transfer model of crown scorch in forest fires. Can. J. For. Res. 36: 2839–2851. doi:10.1139/X06-158.
#' 
#'   Michaletz ST, Johnson EA. 2008. A biophysical process model of tree mortality in surface fires. Canadian Journal of Forest Research 38: 2013–2029.
#'   
#' @name fire_severity
fire_plumeTemperature <- function(Ib_surf, z, T_air = 25.0, rho_air = 1.169) {
    .Call(`_medfate_plumeTemperature`, Ib_surf, z, T_air, rho_air)
}

#' @rdname fire_severity
fire_barkThermalDiffusivity <- function(fmc_bark, rho_bark = 500.0, T_air = 25.0) {
    .Call(`_medfate_barkThermalDiffusivity`, fmc_bark, rho_bark, T_air)
}

#' @rdname fire_severity
fire_radialBoleNecrosis <- function(Ib_surf, t_res, bark_diffusivity, T_air = 25.0, rho_air = 1.169, T_necrosis = 60.0) {
    .Call(`_medfate_radialBoleNecrosis`, Ib_surf, t_res, bark_diffusivity, T_air, rho_air, T_necrosis)
}

#' @rdname fire_severity
fire_leafThermalFactor <- function(SLA, h = 130.0, c = 2500.0) {
    .Call(`_medfate_leafThermalFactor`, SLA, h, c)
}

#' @rdname fire_severity
fire_necrosisCriticalTemperature <- function(t_res, thermal_factor, T_air = 25.0, T_necrosis = 60.0) {
    .Call(`_medfate_necrosisCriticalTemperature`, t_res, thermal_factor, T_air, T_necrosis)
}

#' @rdname fire_severity
fire_necrosisHeight <- function(Ib_surf, t_res, thermal_factor, T_air = 25.0, rho_air = 1.169, T_necrosis = 60.0) {
    .Call(`_medfate_necrosisHeight`, Ib_surf, t_res, thermal_factor, T_air, rho_air, T_necrosis)
}

.treeBasalArea <- function(N, dbh) {
    .Call(`_medfate_treeBasalArea`, N, dbh)
}

.shrubCrownRatio <- function(SP, SpParams) {
    .Call(`_medfate_shrubCrownRatioAllometric`, SP, SpParams)
}

.shrubCover <- function(x, excludeMinHeight = 0.0) {
    .Call(`_medfate_shrubCover`, x, excludeMinHeight)
}

.shrubPhytovolume <- function(SP, Cover, H, SpParams) {
    .Call(`_medfate_shrubPhytovolumeAllometric`, SP, Cover, H, SpParams)
}

#' Woody plant cohort description functions
#'
#' Functions to calculate attributes of woody plants in a \code{\link{forest}} object.
#' 
#' @param x An object of class \code{\link{forest}}.
#' @param SpParams A data frame with species parameters (see \code{\link{SpParamsMED}}).
#' @param parName A string with a parameter name.
#' @param gdd Growth degree days (to account for leaf phenology effects).
#' @param AET Actual annual evapotranspiration (in mm).
#' @param smallBranchDecompositionRate Decomposition rate of small branches.
#' @param includeDead A flag to indicate that standing dead fuels (dead branches) are included.
#' @param treeOffset,shrubOffset Integers to offset cohort IDs.
#' @param fillMissing A boolean flag to try imputation on missing values.
#' @param self_proportion Proportion of the target cohort included in the assessment
#' @param bounded A boolean flag to indicate that extreme values should be prevented (maximum tree LAI = 7 and maximum shrub LAI = 3)
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @return
#' A vector with values for each woody plant cohort of the input \code{\link{forest}} object:
#' \itemize{
#'   \item{\code{plant_basalArea}: Tree basal area (m2/ha).}
#'   \item{\code{plant_largerTreeBasalArea}: Basal area (m2/ha) of trees larger (in diameter) than the tree. Half of the trees of the same record are included.}
#'   \item{\code{plant_characterParameter}: The parameter values of each plant, as strings.}
#'   \item{\code{plant_cover}: Shrub cover (in percent).}
#'   \item{\code{plant_crownBaseHeight}: The height corresponding to the start of the crown (in cm).}
#'   \item{\code{plant_crownLength}: The difference between crown base height and total height (in cm).}
#'   \item{\code{plant_crownRatio}: The ratio between crown length and total height (between 0 and 1).}
#'   \item{\code{plant_density}: Plant density (ind/ha). Tree density is directly taken from the forest object, while the shrub density is estimated from cover and height by calculating the area of a single individual.}
#'   \item{\code{plant_equilibriumLeafLitter}: Litter biomass of leaves at equilibrium (in kg/m2).}
#'   \item{\code{plant_equilibriumSmallBranchLitter}: Litter biomass of small branches (< 6.35 mm diameter) at equilibrium (in kg/m2).}
#'   \item{\code{plant_foliarBiomass}: Standing biomass of leaves (in kg/m2).}
#'   \item{\code{plant_fuelLoading}: Fine fuel load (in kg/m2).}
#'   \item{\code{plant_height}: Total height (in cm).}
#'   \item{\code{plant_ID}: Cohort coding for simulation functions (concatenation of 'T' (Trees) or 'S' (Shrub), cohort index and species index).}
#'   \item{\code{plant_LAI}: Leaf area index (m2/m2).}
#'   \item{\code{plant_individualArea}: Area (m2) occupied by a shrub individual.}
#'   \item{\code{plant_parameter}: The parameter values of each plant, as numeric.}
#'   \item{\code{plant_phytovolume}: Shrub phytovolume (m3/m2).}
#'   \item{\code{plant_species}: Species identity integer (indices start with 0).}
#'   \item{\code{plant_speciesName}: String with species taxonomic name (or a functional group).}
#' }
#' 
#' @seealso  \code{\link{spwb}}, \code{\link{forest}}, \code{\link{summary.forest}}
#'
#' @examples
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Load example plot
#' data(exampleforestMED)
#' 
#' #A plant-level way to obtain stand basal area
#' sum(plant_basalArea(exampleforestMED, SpParamsMED), na.rm=TRUE)
#' 
#' #The analogous plant-level function for LAI
#' sum(plant_LAI(exampleforestMED, SpParamsMED))
#'   
#' #The analogous plant-level function for fuel loading
#' sum(plant_fuelLoading(exampleforestMED, SpParamsMED))
#'       
#' #Summary function for 'forest' objects can be also used
#' summary(exampleforestMED, SpParamsMED)
#' 
#' #Cohort IDs in the models
#' plant_ID(exampleforestMED, SpParamsMED)
#'       
#' @name plant_values
plant_ID <- function(x, SpParams, treeOffset = 0L, shrubOffset = 0L) {
    .Call(`_medfate_cohortIDs`, x, SpParams, treeOffset, shrubOffset)
}

#' @rdname plant_values
plant_basalArea <- function(x, SpParams) {
    .Call(`_medfate_cohortBasalArea`, x, SpParams)
}

#' @rdname plant_values
plant_largerTreeBasalArea <- function(x, SpParams, self_proportion = 0.5) {
    .Call(`_medfate_cohortLargerTreeBasalArea`, x, SpParams, self_proportion)
}

#' @rdname plant_values
plant_cover <- function(x, SpParams) {
    .Call(`_medfate_cohortCover`, x, SpParams)
}

#' @rdname plant_values
plant_species <- function(x, SpParams) {
    .Call(`_medfate_cohortSpecies`, x, SpParams)
}

#' @rdname plant_values
plant_speciesName <- function(x, SpParams) {
    .Call(`_medfate_cohortSpeciesName`, x, SpParams)
}

#' @rdname plant_values
plant_density <- function(x, SpParams) {
    .Call(`_medfate_cohortDensity`, x, SpParams)
}

#' @rdname plant_values
plant_height <- function(x, SpParams) {
    .Call(`_medfate_cohortHeight`, x, SpParams)
}

#' @rdname plant_values
plant_individualArea <- function(x, SpParams) {
    .Call(`_medfate_cohortIndividualArea`, x, SpParams)
}

#' @rdname plant_values
plant_crownRatio <- function(x, SpParams) {
    .Call(`_medfate_cohortCrownRatio`, x, SpParams)
}

#' @rdname plant_values
plant_crownBaseHeight <- function(x, SpParams) {
    .Call(`_medfate_cohortCrownBaseHeight`, x, SpParams)
}

#' @rdname plant_values
plant_crownLength <- function(x, SpParams) {
    .Call(`_medfate_cohortCrownLength`, x, SpParams)
}

#' @rdname plant_values
plant_foliarBiomass <- function(x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_cohortFoliarBiomass`, x, SpParams, gdd)
}

#' @rdname plant_values
plant_fuelLoading <- function(x, SpParams, gdd = NA_real_, includeDead = TRUE) {
    .Call(`_medfate_cohortFuelLoading`, x, SpParams, gdd, includeDead)
}

#' @rdname plant_values
plant_equilibriumLeafLitter <- function(x, SpParams, AET = 800) {
    .Call(`_medfate_cohortEquilibriumLeafLitter`, x, SpParams, AET)
}

#' @rdname plant_values
plant_equilibriumSmallBranchLitter <- function(x, SpParams, smallBranchDecompositionRate = 0.81) {
    .Call(`_medfate_cohortEquilibriumSmallBranchLitter`, x, SpParams, smallBranchDecompositionRate)
}

#' @rdname plant_values
plant_phytovolume <- function(x, SpParams) {
    .Call(`_medfate_cohortPhytovolume`, x, SpParams)
}

#' @rdname plant_values
plant_LAI <- function(x, SpParams, gdd = NA_real_, bounded = TRUE) {
    .Call(`_medfate_cohortLAI`, x, SpParams, gdd, bounded)
}

#' Herbaceous description functions
#'
#' Functions to calculate attributes of the herbaceous component of a \code{\link{forest}} object 
#' @param SpParams A data frame with species parameters (see \code{\link{SpParamsMED}}).
#'  
#' @param x An object of class \code{\link{forest}}.
#' 
#' @return
#' A single scalar:
#' \itemize{
#'   \item{\code{herb_foliarBiomass}: Herbaceous biomass of leaves (in kg/m2).}
#'   \item{\code{herb_fuelLoading}: Herbaceous fine fuel loading (in kg/m2).}
#'   \item{\code{herb_LAI}: Herbaceous leaf area index (m2/m2).}
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{spwb}}, \code{\link{forest}}, \code{\link{plant_basalArea}}, \code{\link{summary.forest}}
#' 
#' @name herb_values
herb_foliarBiomass <- function(x, SpParams) {
    .Call(`_medfate_herbFoliarBiomass`, x, SpParams)
}

#' @rdname herb_values
herb_fuelLoading <- function(x, SpParams) {
    .Call(`_medfate_herbFuelLoading`, x, SpParams)
}

#' @rdname herb_values
herb_LAI <- function(x, SpParams) {
    .Call(`_medfate_herbLAI`, x, SpParams)
}

#' Species description functions
#'
#' Functions to calculate attributes of a \code{\link{forest}} object by species or to extract species parameters from a species parameter table (\code{\link{SpParamsMED}}).
#' 
#' @param x An object of class \code{\link{forest}}.
#' @param SpParams A data frame with species parameters (see \code{\link{SpParamsMED}}).
#' @param gdd Growth degree days (to account for leaf phenology effects).
#' @param includeDead A flag to indicate that standing dead fuels (dead branches) are included.
#' @param species A character vector of species names.
#' @param parName A string with a parameter name.
#' @param fillMissing A boolean flag to try imputation on missing values.
#' @param bounded A boolean flag to indicate that extreme values should be prevented (maximum tree LAI = 7 and maximum shrub LAI = 3)
#' 
#' @return
#' A vector with values for each species in \code{SpParams}:
#' \itemize{
#'   \item{\code{species_basalArea}: Species basal area (m2/ha).}
#'   \item{\code{species_cover}: Shrub cover (in percent).}
#'   \item{\code{species_density}: Plant density (ind/ha). Tree density is directly taken from the forest object, while the shrub density is estimated from cover and height by calculating the area of a single individual.}
#'   \item{\code{species_foliarBiomass}: Standing biomass of leaves (in kg/m2).}
#'   \item{\code{species_fuel}: Fine fuel load (in kg/m2).}
#'   \item{\code{species_LAI}: Leaf area index (m2/m2).}
#'   \item{\code{species_phytovolume}: Shrub phytovolume (m3/m2).}
#'   \item{\code{species_parameter}: A numeric vector with the parameter values of each input species.}
#'   \item{\code{species_characterParameter}: A character vector with the parameter values of each input species.}
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{spwb}}, \code{\link{forest}}, \code{\link{plant_basalArea}}, \code{\link{summary.forest}}
#' 
#' @examples
#' # Default species parameterization
#' data(SpParamsMED)
#' 
#' # Load example plot
#' data(exampleforestMED)
#' 
#' # Species basal area in the forest plot
#' species_basalArea(exampleforestMED, SpParamsMED)
#'   
#' # Value of parameter "Psi_Extract" for two species
#' species_parameter(c("Pinus halepensis", "Quercus ilex"), SpParamsMED, "Psi_Extract")
#'     
#' @name species_values
species_basalArea <- function(x, SpParams) {
    .Call(`_medfate_speciesBasalArea`, x, SpParams)
}

#' @rdname species_values
species_cover <- function(x, SpParams) {
    .Call(`_medfate_speciesCover`, x, SpParams)
}

#' @rdname species_values
species_density <- function(x, SpParams) {
    .Call(`_medfate_speciesDensity`, x, SpParams)
}

#' @rdname species_values
species_foliarBiomass <- function(x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_speciesFoliarBiomass`, x, SpParams, gdd)
}

#' @rdname species_values
species_fuelLoading <- function(x, SpParams, gdd = NA_real_, includeDead = TRUE) {
    .Call(`_medfate_speciesFuelLoading`, x, SpParams, gdd, includeDead)
}

#' @rdname species_values
species_LAI <- function(x, SpParams, gdd = NA_real_, bounded = TRUE) {
    .Call(`_medfate_speciesLAI`, x, SpParams, gdd, bounded)
}

#' @rdname stand_values
stand_basalArea <- function(x, minDBH = 7.5) {
    .Call(`_medfate_standBasalArea`, x, minDBH)
}

#' @rdname stand_values
stand_foliarBiomass <- function(x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_standFoliarBiomass`, x, SpParams, gdd)
}

#' @rdname stand_values
stand_fuelLoading <- function(x, SpParams, gdd = NA_real_, includeDead = TRUE) {
    .Call(`_medfate_standFuelLoading`, x, SpParams, gdd, includeDead)
}

#' @rdname stand_values
stand_LAI <- function(x, SpParams, gdd = NA_real_, bounded = TRUE) {
    .Call(`_medfate_standLAI`, x, SpParams, gdd, bounded)
}

.LAIdistributionVectors <- function(z, LAI, H, CR) {
    .Call(`_medfate_LAIdistributionVectors`, z, LAI, H, CR)
}

.LAIdistribution <- function(z, x, SpParams, gdd = NA_real_, bounded = TRUE) {
    .Call(`_medfate_LAIdistribution`, z, x, SpParams, gdd, bounded)
}

.LAIprofileVectors <- function(z, LAI, H, CR) {
    .Call(`_medfate_LAIprofileVectors`, z, LAI, H, CR)
}

.LAIprofile <- function(z, x, SpParams, gdd = NA_real_, bounded = TRUE) {
    .Call(`_medfate_LAIprofile`, z, x, SpParams, gdd, bounded)
}

#' @rdname modelInput
forest2aboveground <- function(x, SpParams, gdd = NA_real_, loading = FALSE) {
    .Call(`_medfate_forest2aboveground`, x, SpParams, gdd, loading)
}

#' @rdname modelInput
forest2belowground <- function(x, soil, SpParams) {
    .Call(`_medfate_forest2belowground`, x, soil, SpParams)
}

.fuelConditions <- function(airTemp, airHumidity, fuelRadiation, fuelWindSpeed) {
    .Call(`_medfate_fuelConditions`, airTemp, airHumidity, fuelRadiation, fuelWindSpeed)
}

.EMCdesorption <- function(fuelTemperature, fuelHumidity) {
    .Call(`_medfate_EMCdesorption`, fuelTemperature, fuelHumidity)
}

.EMCadsorption <- function(fuelTemperature, fuelHumidity) {
    .Call(`_medfate_EMCadsorption`, fuelTemperature, fuelHumidity)
}

.EMCSimard <- function(fuelTemperature, fuelHumidity) {
    .Call(`_medfate_EMCSimard`, fuelTemperature, fuelHumidity)
}

.woodyFuelProfile <- function(z, x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_woodyFuelProfile`, z, x, SpParams, gdd)
}

.layerCohortFuelLoading <- function(minHeight, maxHeight, cohortLoading, H, CR) {
    .Call(`_medfate_layerCohortFuelLoading`, minHeight, maxHeight, cohortLoading, H, CR)
}

.layerFuelLoading <- function(minHeight, maxHeight, cohortLoading, H, CR) {
    .Call(`_medfate_layerFuelLoading`, minHeight, maxHeight, cohortLoading, H, CR)
}

.layerLAI <- function(minHeight, maxHeight, cohortLAI, H, CR) {
    .Call(`_medfate_layerLAI`, minHeight, maxHeight, cohortLAI, H, CR)
}

.layerFuelAverageSpeciesParameter <- function(spParName, minHeight, maxHeight, x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_layerFuelAverageSpeciesParameter`, spParName, minHeight, maxHeight, x, SpParams, gdd)
}

.layerFuelAverageParameter <- function(minHeight, maxHeight, cohortParameter, cohortLoading, H, CR) {
    .Call(`_medfate_layerFuelAverageParameter`, minHeight, maxHeight, cohortParameter, cohortLoading, H, CR)
}

.layerFuelAverageCrownLength <- function(minHeight, maxHeight, cohortCrownLength, cohortLoading, H, CR) {
    .Call(`_medfate_layerFuelAverageCrownLength`, minHeight, maxHeight, cohortCrownLength, cohortLoading, H, CR)
}

#' Fuel stratification and fuel characteristics
#' 
#' Function \code{fuel_stratification} provides a stratification of the stand into understory and canopy strata. 
#' Function \code{fuel_FCCS} calculates fuel characteristics from a \code{forest} object 
#' following an adaptation of the protocols described for the Fuel Characteristics Classification System (Prichard et al. 2013). 
#' Function \code{fuel_windAdjustmentFactor} determines the adjustment factor of wind for surface fires, according to Andrews (2012). 
#' 
#' @param object An object of class \code{\link{forest}}
#' @param SpParams A data frame with species parameters (see \code{\link{SpParamsMED}}).
#' @param cohortFMC A numeric vector of (actual) fuel moisture content by cohort.
#' @param gdd Growth degree-days.
#' @param heightProfileStep Precision for the fuel bulk density profile.
#' @param maxHeightProfile Maximum height for the fuel bulk density profile.
#' 
#' @return 
#' Function \code{fuel_FCCS} returns a data frame with five rows corresponding to  fuel layers: 
#' \code{canopy}, \code{shrub}, \code{herb}, \code{woody} and \code{litter}. Columns correspond fuel properties:
#' \itemize{
#'   \item{\code{w}: Fine fuel loading (in kg/m2).}
#'   \item{\code{cover}: Percent cover.}
#'   \item{\code{hbc}: Height to base of crowns (in m).}
#'   \item{\code{htc}: Height to top of crowns (in m).}
#'   \item{\code{delta}: Fuel depth (in m).}
#'   \item{\code{rhob}: Fuel bulk density (in kg/m3).}
#'   \item{\code{rhop}: Fuel particle density (in kg/m3).}
#'   \item{\code{PV}: Particle volume (in m3/m2).}
#'   \item{\code{beta}: Packing ratio (unitless).}
#'   \item{\code{betarel}: Relative packing ratio (unitless).}
#'   \item{\code{etabetarel}: Reaction efficiency (unitless).}
#'   \item{\code{sigma}: Surface area-to-volume ratio (m2/m3).}
#'   \item{\code{pDead}: Proportion of dead fuels.}
#'   \item{\code{FAI}: Fuel area index (unitless).}
#'   \item{\code{h}: High heat content (in kJ/kg).}
#'   \item{\code{RV}: Reactive volume (in m3/m2).}
#'   \item{\code{MinFMC}: Minimum fuel moisture content (as percent over dry weight).}
#'   \item{\code{MaxFMC}: Maximum fuel moisture content (as percent over dry weight).}
#'   \item{\code{ActFMC}: Actual fuel moisture content (as percent over dry weight). These are set to \code{NA} if parameter \code{cohortFMC} is empty.}
#' }
#' 
#' Function \code{fuel_stratification} returns a list with the following items:
#'   \itemize{
#'     \item{\code{surfaceLayerBaseHeight}: Base height of crowns of shrubs in the surface layer (in cm).}
#'     \item{\code{surfaceLayerTopHeight}: Top height of crowns of shrubs in the surface layer (in cm).}
#'     \item{\code{understoryLAI}: Cumulated LAI of the understory layer (i.e. leaf area comprised between surface layer base and top heights).}
#'     \item{\code{canopyBaseHeight}: Base height of tree crowns in the canopy (in cm).}
#'     \item{\code{canopyTopHeight}: Top height of tree crowns in the canopy (in cm).}
#'     \item{\code{canopyLAI}: Cumulated LAI of the canopy (i.e. leaf area comprised between canopy base and top heights).}
#'   }
#'   
#' Function \code{fuel_cohortFineFMC} returns a list with three matrices (for leaves, twigs and fine fuels). 
#' Each of them contains live moisture content values for each day (in rows) and plant cohort (in columns).
#' 
#' Function \code{fuel_windAdjustmentFactor} returns a value between 0 and 1.
#' 
#' @references
#' Andrews, P. L. 2012. Modeling wind adjustment factor and midflame wind speed for Rothermel’s surface fire spread model. USDA Forest Service - General Technical Report RMRS-GTR:1–39.
#' 
#' Prichard, S. J., D. V Sandberg, R. D. Ottmar, E. Eberhardt, A. Andreu, P. Eagle, and K. Swedin. 2013. Classification System Version 3.0: Technical Documentation.
#' 
#' Reinhardt, E., D. Lutes, and J. Scott. 2006. FuelCalc: A method for estimating fuel characteristics. Pages 273–282.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{fire_FCCS}}, \code{\link{spwb}}
#' 
#' @examples
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Show stratification of fuels
#' fuel_stratification(exampleforestMED, SpParamsMED)
#'   
#' #Calculate fuel properties according to FCCS
#' fccs = fuel_FCCS(exampleforestMED, SpParamsMED)
#' fccs
#' 
#' fuel_windAdjustmentFactor(fccs$htc[2], fccs$hbc[1], fccs$htc[1], fccs$cover[1])
#' 
#' @name fuel_properties
fuel_stratification <- function(object, SpParams, gdd = NA_real_, heightProfileStep = 10.0, maxHeightProfile = 5000.0, bulkDensityThreshold = 0.05) {
    .Call(`_medfate_fuelLiveStratification`, object, SpParams, gdd, heightProfileStep, maxHeightProfile, bulkDensityThreshold)
}

#' @rdname fuel_properties
#' 
#' @param bulkDensityThreshold Minimum fuel bulk density to delimit fuel strata.
#' @param depthMode Specifies how fuel depth (and therefore canopy and understory bulk density) should be estimated: 
#'   \itemize{
#'     \item{\code{"crownaverage"}: As weighed average of crown lengths using loadings as weights.}
#'     \item{\code{"profile"}: As the difference of base and top heights in bulk density profiles.}  
#'     \item{\code{"absoluteprofile"}: As the difference of absolute base and absolute top heights in bulk density profiles.}  
#'   }
fuel_FCCS <- function(object, SpParams, cohortFMC = as.numeric( c()), gdd = NA_real_, heightProfileStep = 10.0, maxHeightProfile = 5000, bulkDensityThreshold = 0.05, depthMode = "crownaverage") {
    .Call(`_medfate_FCCSproperties`, object, SpParams, cohortFMC, gdd, heightProfileStep, maxHeightProfile, bulkDensityThreshold, depthMode)
}

#' Mortality
#' 
#' A simple sigmoid function to determine a daily mortality likelihood according 
#' to the value of a stress variable.
#'
#' @param stressValue Current value of the stress variable (0 to 1, 
#'                    with higher values indicate stronger stress).
#' @param stressThreshold Threshold to indicate 50\% annual mortality probability.
#' 
#' @return Returns a probability (between 0 and 1)
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{growth}}, \code{\link{regeneration}}
#' 
mortality_dailyProbability <- function(stressValue, stressThreshold) {
    .Call(`_medfate_dailyMortalityProbability`, stressValue, stressThreshold)
}

#' @rdname spwb_day
growth_day <- function(x, date, meteovec, latitude, elevation, slope, aspect, runon = 0.0, modifyInput = TRUE) {
    .Call(`_medfate_growthDay`, x, date, meteovec, latitude, elevation, slope, aspect, runon, modifyInput)
}

#' Forest growth
#' 
#' Function \code{growth} is a process-based model that performs energy, water and carbon balances; 
#' and determines changes in water/carbon pools, functional variables (leaf area, sapwood area, root area) 
#' and structural ones (tree diameter, tree height, shrub cover) for woody plant cohorts in a given forest stand 
#' during a period specified in the input climatic data. 
#' 
#' @param x An object of class \code{\link{growthInput}}.
#' @param meteo A data frame with daily meteorological data series (see \code{\link{spwb}}).
#' @param latitude Latitude (in degrees).
#' @param elevation,slope,aspect Elevation above sea level (in m), slope (in degrees) and aspect (in degrees from North). 
#' @param CO2ByYear A named numeric vector with years as names and atmospheric CO2 concentration (in ppm) as values. Used to specify annual changes in CO2 concentration along the simulation (as an alternative to specifying daily values in \code{meteo}).
#' 
#' @details
#' Detailed model description is available in the medfate book. 
#' Simulations using the 'Sperry' or 'Cochard' transpiration modes are computationally much more expensive 
#' than those using the 'Granier' transpiration mode. 
#' 
#' @return
#' A list of class 'growth' with the following elements:
#' \itemize{
#'   \item{\code{"latitude"}: Latitude (in degrees) given as input.} 
#'   \item{\code{"topography"}: Vector with elevation, slope and aspect given as input.} 
#'   \item{\code{"weather"}: A copy of the input weather data frame.}
#'   \item{\code{"growthInput"}: A copy of the object \code{x} of class \code{\link{growthInput}} given as input.}
#'   \item{\code{"growthOutput"}: An copy of the final state of the object \code{x} of class \code{\link{growthInput}}.}
#'   \item{\code{"WaterBalance"}: A data frame where different water balance variables (see \code{\link{spwb}}).}
#'   \item{\code{"EnergyBalance"}: A data frame with the daily values of energy balance components for the soil and the canopy (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}; see \code{\link{spwb}}).}
#'   \item{\code{"CarbonBalance"}: A data frame where different stand-level carbon balance components (gross primary production, maintenance respiration, synthesis respiration and net primary production), all in g C · m-2.}
#'   \item{\code{"BiomassBalance"}: A data frame with the daily values of stand biomass balance components (in g dry · m-2.}
#'   \item{\code{"Temperature"}: A data frame with the daily values of minimum/mean/maximum temperatures for the atmosphere (input), canopy and soil (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}; see \code{\link{spwb}}).}
#'   \item{\code{"Soil"}: A data frame where different soil variables  (see \code{\link{spwb}}).}
#'   \item{\code{"Stand"}: A data frame where different stand-level variables (see \code{\link{spwb}}).}
#'   \item{\code{"Plants"}: A list of daily results for plant cohorts (see \code{\link{spwb}}).}
#'   \item{\code{"SunlitLeaves"} and \code{"ShadeLeaves"}: A list with daily results for sunlit and shade leaves (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}; see \code{\link{spwb}}).}
#'   \item{\code{"LabileCarbonBalance"}: A list of daily labile carbon balance results for plant cohorts, with elements:}
#'   \itemize{
#'     \item{\code{"GrossPhotosynthesis"}: Daily gross photosynthesis per dry weight of living biomass (g gluc · g dry-1).}
#'     \item{\code{"MaintentanceRespiration"}: Daily maintenance respiration per dry weight of living biomass (g gluc · g dry-1).}
#'     \item{\code{"GrowthCosts"}: Daily growth costs per dry weight of living biomass (g gluc · g dry-1).}
#'     \item{\code{"RootExudation"}: Root exudation per dry weight of living biomass (g gluc · g dry-1).}    
#'     \item{\code{"LabileCarbonBalance"}: Daily labile carbon balance (photosynthesis - maintenance respiration - growth costs - root exudation) per dry weight of living biomass (g gluc · g dry-1).}
#'     \item{\code{"SugarLeaf"}: Sugar concentration (mol·l-1) in leaves.}
#'     \item{\code{"StarchLeaf"}: Starch concentration (mol·l-1) in leaves.}
#'     \item{\code{"SugarSapwood"}: Sugar concentration (mol·l-1) in sapwood.}
#'     \item{\code{"StarchSapwood"}: Starch concentration (mol·l-1) in sapwood.}
#'     \item{\code{"SugarTransport"}:  Average instantaneous rate of carbon transferred between leaves and stem compartments via floem (mol gluc·s-1).}
#'   }
#'   \item{\code{"PlantBiomassBalance"}: A list of daily plant biomass balance results for plant cohorts, with elements:}
#'   \itemize{
#'     \item{\code{"StructuralBiomassBalance"}: Daily structural biomass balance (g dry · m-2).}
#'     \item{\code{"LabileBiomassBalance"}: Daily labile biomass balance (g dry · m-2).}
#'     \item{\code{"PlantBiomassBalance"}: Daily plant biomass balance, i.e. labile change + structural change (g dry · m-2).}
#'     \item{\code{"MortalityBiomassLoss"}: Biomass loss due to mortality (g dry · m-2).}    
#'     \item{\code{"CohortBiomassBalance"}: Daily cohort biomass balance (including mortality) (g dry · m-2).}
#'   }
#'   \item{\code{"PlantStructure"}: A list of daily area and biomass values for compartments of plant cohorts, with elements:}
#'   \itemize{
#'     \item{\code{"LeafBiomass"}: Daily amount of leaf structural biomass (in g dry) for an average individual of each plant cohort.}
#'     \item{\code{"SapwoodBiomass"}: Daily amount of sapwood structural biomass (in g dry) for an average individual of each plant cohort.}
#'     \item{\code{"FineRootBiomass"}: Daily amount of fine root biomass (in g dry) for an average individual of each plant cohort.}
#'     \item{\code{"LeafArea"}: Daily amount of leaf area (in m2) for an average individual of each plant cohort.}
#'     \item{\code{"SapwoodArea"}: Daily amount of sapwood area (in cm2) for an average individual of each plant cohort.}
#'     \item{\code{"FineRootArea"}: Daily amount of fine root area (in m2) for an average individual of each plant cohort.}
#'     \item{\code{"HuberValue"}: The ratio of sapwood area to (target) leaf area (in cm2/m2).}
#'     \item{\code{"RootAreaLeafArea"}: The ratio of fine root area to (target) leaf area (in m2/m2).}
#'     \item{\code{"DBH"}: Diameter at breast height (in cm) for an average individual of each plant cohort.}
#'     \item{\code{"Height"}: Height (in cm) for an average individual of each plant cohort.}
#'   }
#'   \item{\code{"GrowthMortality"}: A list of daily growth and mortality rates for plant cohorts, with elements:}
#'   \itemize{
#'     \item{\code{"LAgrowth"}: Leaf area growth (in m2·day-1) for an average individual of each plant cohort.}
#'     \item{\code{"SAgrowth"}: Sapwood area growth rate (in cm2·day-1) for an average individual of each plant cohort.}
#'     \item{\code{"FRAgrowth"}: Fine root area growth (in m2·day-1) for an average individual of each plant cohort.}
#'     \item{\code{"StarvationRate"}: Daily mortality rate from starvation (ind/d-1).}
#'     \item{\code{"DessicationRate"}: Daily mortality rate from dessication (ind/d-1).}
#'     \item{\code{"MortalityRate"}: Daily mortality rate (any cause) (ind/d-1).}
#'   }
#'   \item{\code{"subdaily"}: A list of objects of class \code{\link{growth_day}}, one per day simulated (only if required in \code{control} parameters, see \code{\link{defaultControl}}).}
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{growthInput}}, \code{\link{growth_day}}, \code{\link{plot.growth}}
#' 
#' @references
#' De Cáceres M, Molowny-Horas R, Cabon A, Martínez-Vilalta J, Mencuccini M, García-Valdés R, Nadal-Sala D, Sabaté S, 
#' Martin-StPaul N, Morin X, D'Adamo F, Batllori E, Améztegui A (2023) MEDFATE 2.9.3: A trait-enabled model to simulate 
#' Mediterranean forest function and dynamics at regional scales. 
#' Geoscientific Model Development 16: 3165-3201 (https://doi.org/10.5194/gmd-16-3165-2023).
#' 
#' @examples
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#'   
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#'   
#' #Initialize soil with default soil params (4 layers)
#' examplesoil <- soil(defaultSoilParams(4))
#' 
#' #Initialize vegetation input
#' x1 <- forest2growthInput(exampleforestMED, examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' G1 <- growth(x1, examplemeteo, latitude = 41.82592, elevation = 100)
#'  
#' \donttest{
#' #Switch to 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' 
#' #Initialize vegetation input
#' x2 <- forest2growthInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' G2 <-growth(x2, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' #Switch to 'Cochard' transpiration mode
#' control <- defaultControl("Cochard")
#' 
#' #Initialize vegetation input
#' x3 <- forest2growthInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' G3 <-growth(x3, examplemeteo, latitude = 41.82592, elevation = 100)
#' }
#'       
growth <- function(x, meteo, latitude, elevation = NA_real_, slope = NA_real_, aspect = NA_real_, CO2ByYear = numeric(0)) {
    .Call(`_medfate_growth`, x, meteo, latitude, elevation, slope, aspect, CO2ByYear)
}

#' Hydraulic confuctance functions
#' 
#' Set of functions used in the calculation of soil and plant hydraulic conductance.
#'
#' @param psi A scalar (or a vector, depending on the function) with water potential (in MPa).
#' @param K Whole-plant relative conductance (0-1).
#' @param psi_extract Soil water potential (in MPa) corresponding to 50\% whole-plant relative transpiration.
#' @param exp_extract Exponent of the whole-plant relative transpiration Weibull function.
#' @param v Proportion of fine roots within each soil layer.
#' @param krhizomax Maximum rhizosphere hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param kxylemmax Maximum xylem hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param c,d Parameters of the Weibull function (generic xylem vulnerability curve).
#' @param n,alpha Parameters of the Van Genuchten function (rhizosphere vulnerability curve).
#' @param kxylem Xylem hydraulic conductance (defined as flow per surface unit and per pressure drop).
#' @param pCrit Proportion of maximum conductance considered critical for hydraulic functioning.
#' @param psi50,psi88,psi12 Water potentials (in MPa) corresponding to 50\%, 88\% and 12\% percent conductance loss.
#' @param temp Temperature (in degrees Celsius).
#' 
#' @details Details of plant hydraulic models are given the medfate book. 
#' Function \code{hydraulics_vulnerabilityCurvePlot} draws a plot of the vulnerability curves for the given \code{soil} object and network properties of each plant cohort in \code{x}.
#' 
#' @return
#' Values returned for each function are:
#' \itemize{
#'   \item{\code{hydraulics_psi2K}: Whole-plant relative conductance (0-1).}
#'   \item{\code{hydraulics_K2Psi}: Soil water potential (in MPa) corresponding to the given whole-plant relative conductance value (inverse of \code{hydraulics_psi2K()}).}
#'   \item{\code{hydraulics_averagePsi}: The average water potential (in MPa) across soil layers.}
#'   \item{\code{hydraulics_vanGenuchtenConductance}: Rhizosphere conductance corresponding to an input water potential (soil vulnerability curve).}
#'   \item{\code{hydraulics_xylemConductance}: Xylem conductance (flow rate per pressure drop) corresponding to an input water potential (plant vulnerability curve).}
#'   \item{\code{hydraulics_xylemPsi}: Xylem water potential (in MPa) corresponding to an input xylem conductance (flow rate per pressure drop).}
#'   \item{\code{hydraulics_psi2Weibull}: Parameters of the Weibull vulnerability curve that goes through the supplied psi50 and psi88 values.}
#' }
#' 
#' @references
#' Sperry, J. S., F. R. Adler, G. S. Campbell, and J. P. Comstock. 1998. Limitation of plant water use by rhizosphere and xylem conductance: results from a model. Plant, Cell and Environment 21:347–359.
#' 
#' Sperry, J. S., and D. M. Love. 2015. What plant hydraulics can tell us about responses to climate-change droughts. New Phytologist 207:14–27.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{hydraulics_supplyFunctionPlot}}, \code{\link{hydraulics_maximumStemHydraulicConductance}}, \code{\link{spwb}}, \code{\link{soil}}
#' 
#' @examples
#' 
#' #Manual display of vulnerability curve
#' kstemmax = 4 # in mmol·m-2·s-1·MPa-1
#' stemc = 3 
#' stemd = -4 # in MPa
#' psiVec = seq(-0.1, -7.0, by =-0.01)
#' kstem = unlist(lapply(psiVec, hydraulics_xylemConductance, kstemmax, stemc, stemd))
#' plot(-psiVec, kstem, type="l",ylab="Xylem conductance (mmol·m-2·s-1·MPa-1)", 
#'      xlab="Canopy pressure (-MPa)", lwd=1.5,ylim=c(0,kstemmax))
#' 
#' #Load example dataset
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Initialize soil with default soil params (2 layers)
#' examplesoil = soil(defaultSoilParams(2)) 
#' 
#' #Initialize control parameters
#' control = defaultControl("Granier")
#' 
#' #Switch to 'Sperry' transpiration mode
#' control = defaultControl("Sperry")
#' 
#' #Initialize input
#' x = forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' #Leaf vulnerability curves
#' hydraulics_vulnerabilityCurvePlot(x, type="leaf")
#' 
#' #Stem vulnerability curves
#' hydraulics_vulnerabilityCurvePlot(x, type="stem")
#'              
#' @name hydraulics_conductancefunctions
hydraulics_psi2K <- function(psi, psi_extract, exp_extract = 3.0) {
    .Call(`_medfate_Psi2K`, psi, psi_extract, exp_extract)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_K2Psi <- function(K, psi_extract, exp_extract = 3.0) {
    .Call(`_medfate_K2Psi`, K, psi_extract, exp_extract)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_averagePsi <- function(psi, v, exp_extract, psi_extract) {
    .Call(`_medfate_averagePsi`, psi, v, exp_extract, psi_extract)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_xylemConductance <- function(psi, kxylemmax, c, d) {
    .Call(`_medfate_xylemConductance`, psi, kxylemmax, c, d)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_xylemPsi <- function(kxylem, kxylemmax, c, d) {
    .Call(`_medfate_xylemPsi`, kxylem, kxylemmax, c, d)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_psiCrit <- function(c, d, pCrit = 0.001) {
    .Call(`_medfate_psiCrit`, c, d, pCrit)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_vanGenuchtenConductance <- function(psi, krhizomax, n, alpha) {
    .Call(`_medfate_vanGenuchtenConductance`, psi, krhizomax, n, alpha)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_correctConductanceForViscosity <- function(kxylem, temp) {
    .Call(`_medfate_correctConductanceForViscosity`, kxylem, temp)
}

#' @rdname hydraulics_conductancefunctions
hydraulics_psi2Weibull <- function(psi50, psi88 = NA_real_, psi12 = NA_real_) {
    .Call(`_medfate_psi2Weibull`, psi50, psi88, psi12)
}

.Egamma <- function(psi, kxylemmax, c, d, psiCav = 0.0) {
    .Call(`_medfate_Egamma`, psi, kxylemmax, c, d, psiCav)
}

.Egammainv <- function(Eg, kxylemmax, c, d, psiCav = 0.0) {
    .Call(`_medfate_Egammainv`, Eg, kxylemmax, c, d, psiCav)
}

#' Hydraulic supply functions
#' 
#' Set of functions used in the implementation of hydraulic supply functions (Sperry and Love 2015).
#'
#' @param v Proportion of fine roots within each soil layer.
#' @param krhizomax Maximum rhizosphere hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param kxylemmax Maximum xylem hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param kleafmax Maximum leaf hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param kstemmax Maximum stem xylem hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param E Flow per surface unit.
#' @param Emax Maximum flow per surface unit.
#' @param Erootcrown Flow per surface unit at the root crown.
#' @param psi Water potential (in MPa).
#' @param psiPrev Water potential (in MPa) in the previous time step.
#' @param psiDownstream Water potential upstream (in MPa).
#' @param psiUpstream Water potential upstream (in MPa). In a one-component model corresponds to soil potential. In a two-component model corresponds to the potential inside the roots.
#' @param psiCav Minimum water potential (in MPa) experienced (for irreversible cavitation).
#' @param minFlow Minimum flow in supply function.
#' @param psiPlant Plant water potential (in MPa).
#' @param hydraulicNetwork List with the hydraulic characteristics of nodes in the hydraulic network.
#' @param psiFineRoot Water potential (in MPa) inside fine roots.
#' @param psiSoil Soil water potential (in MPa). A scalar or a vector depending on the function.
#' @param psiRhizo Soil water potential (in MPa) in the rhizosphere (root surface).
#' @param psiRootCrown Soil water potential (in MPa) at the root crown.
#' @param psiStep Water potential precision (in MPa).
#' @param psiTol Precision for water potential estimates (in MPa).
#' @param psiIni Vector of initial water potential values (in MPa).
#' @param psiMax Minimum (maximum in absolute value) water potential to be considered (in MPa).
#' @param pCrit Critical water potential (in MPa).
#' @param PLCprev Previous proportion of loss conductance [0-1].
#' @param V Capacity of the compartment per leaf area (in L/m2).
#' @param fapo Apoplastic fraction (proportion) in the segment.
#' @param pi0 Full turgor osmotic potential (MPa).
#' @param eps Bulk modulus of elasticity (MPa).
#' @param dE Increment of flow per surface unit.
#' @param ETol Precision for water flow per surface unit.
#' @param c,d Parameters of the Weibull function (generic xylem vulnerability curve).
#' @param stemc,stemd Parameters of the Weibull function for stems (stem xylem vulnerability curve).
#' @param leafc,leafd Parameters of the Weibull function for leaves (leaf vulnerability curve).
#' @param n,alpha,l Parameters of the Van Genuchten function (rhizosphere vulnerability curve).
#' @param allowNegativeFlux A boolean to indicate wether negative flux (i.e. from plant to soil) is allowed.
#' @param maxNsteps Maximum number of steps in the construction of supply functions.
#' @param ntrial Maximum number of steps in Newton-Raphson optimization.
#' @param timestep Time step in seconds.
#' 
#' @details 
#' Function \code{hydraulics_supplyFunctionPlot} draws a plot of the supply function for the given \code{soil} object and network properties of each plant cohort in \code{x}. Function \code{hydraulics_vulnerabilityCurvePlot} draws a plot of the vulnerability curves for the given \code{soil} object and network properties of each plant cohort in \code{x}.
#' 
#' @return
#' Values returned for each function are:
#' \itemize{
#'   \item{\code{hydraulics_E2psiXylem}: The plant (leaf) water potential (in MPa) corresponding to the input flow, according to the xylem supply function and given an upstream (soil or root) water potential.}
#'   \item{\code{hydraulics_E2psiVanGenuchten}: The root water potential (in MPa) corresponding to the input flow, according to the rhizosphere supply function and given a soil water potential.}
#'   \item{\code{hydraulics_E2psiTwoElements}: The plant (leaf) water potential (in MPa) corresponding to the input flow, according to the rhizosphere and plant supply functions and given an input soil water potential.}
#'   \item{\code{hydraulics_E2psiNetwork}: The rhizosphere, root crown and plant (leaf water potential (in MPa) corresponding to the input flow, according to the vulnerability curves of rhizosphere, root and stem elements in a network.}
#'   \item{\code{hydraulics_Ecrit}: The critical flow according to the xylem supply function and given an input soil water potential.}
#'   \item{\code{hydraulics_EVanGenuchten}: The flow (integral of the vulnerability curve) according to the rhizosphere supply function and given an input drop in water potential (soil and rhizosphere).}
#'   \item{\code{hydraulics_EXylem}: The flow (integral of the vulnerability curve) according to the xylem supply function and given an input drop in water potential (rhizosphere and plant).}
#'   \item{\code{hydraulics_supplyFunctionOneXylem}, \code{hydraulics_supplyFunctionTwoElements} and
#'     \code{hydraulics_supplyFunctionNetwork}: A list with different numeric vectors with information of the two-element supply function:
#'     \itemize{
#'       \item{\code{E}: Flow values (supply values).}
#'       \item{\code{FittedE}: Fitted flow values (for \code{hydraulics_supplyFunctionTwoElements}).}
#'       \item{\code{Elayers}: Flow values across the roots of each soil layer (only for \code{hydraulics_supplyFunctionNetwork}).}
#'       \item{\code{PsiRhizo}: Water potential values at the root surface (only for \code{hydraulics_supplyFunctionNetwork}).}
#'       \item{\code{PsiRoot}: Water potential values inside the root crown (not for \code{hydraulics_supplyFunctionOneXylem}).}
#'       \item{\code{PsiPlant}: Water potential values at the canopy (leaf).}
#'       \item{\code{dEdP}: Derivatives of the supply function.}
#'     }
#'   }
#'   \item{\code{hydraulics_supplyFunctionPlot}: If \code{draw = FALSE} a list with the result of calling \code{hydraulics_supplyFunctionNetwork} for each cohort. }
#'   \item{\code{hydraulics_regulatedPsiXylem}: Plant water potential after regulation (one-element loss function) given an input water potential.}
#'   \item{\code{hydraulics_regulatedPsiTwoElements}: Plant water potential after regulation (two-element loss function) given an input soil water potential.}
#' }
#' 
#' @references
#' Sperry, J. S., F. R. Adler, G. S. Campbell, and J. P. Comstock. 1998. Limitation of plant water use by rhizosphere and xylem conductance: results from a model. Plant, Cell and Environment 21:347–359.
#' 
#' Sperry, J. S., and D. M. Love. 2015. What plant hydraulics can tell us about responses to climate-change droughts. New Phytologist 207:14–27.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{hydraulics_psi2K}}, \code{\link{hydraulics_maximumStemHydraulicConductance}}, \code{\link{spwb}}, \code{\link{soil}}
#' 
#' @examples
#' kstemmax = 4 # in mmol·m-2·s-1·MPa-1
#' stemc = 3 
#' stemd = -4 # in MPa
#' psiVec = seq(-0.1, -7.0, by =-0.01)
#' 
#' #Vulnerability curve
#' kstem = unlist(lapply(psiVec, hydraulics_xylemConductance, kstemmax, stemc, stemd))
#' plot(-psiVec, kstem, type="l",ylab="Xylem conductance (mmol·m-2·s-1·MPa-1)", 
#'      xlab="Canopy pressure (-MPa)", lwd=1.5,ylim=c(0,kstemmax))
#' 
#' @name hydraulics_supplyfunctions
hydraulics_EXylem <- function(psiPlant, psiUpstream, kxylemmax, c, d, allowNegativeFlux = TRUE, psiCav = 0.0) {
    .Call(`_medfate_EXylem`, psiPlant, psiUpstream, kxylemmax, c, d, allowNegativeFlux, psiCav)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiXylem <- function(E, psiUpstream, kxylemmax, c, d, psiCav = 0.0) {
    .Call(`_medfate_E2psiXylem`, E, psiUpstream, kxylemmax, c, d, psiCav)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiXylemUp <- function(E, psiDownstream, kxylemmax, c, d, psiCav = 0.0) {
    .Call(`_medfate_E2psiXylemUp`, E, psiDownstream, kxylemmax, c, d, psiCav)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_EVanGenuchten <- function(psiRhizo, psiSoil, krhizomax, n, alpha, l = 0.5) {
    .Call(`_medfate_EVanGenuchten`, psiRhizo, psiSoil, krhizomax, n, alpha, l)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_ECrit <- function(psiUpstream, kxylemmax, c, d, pCrit = 0.001) {
    .Call(`_medfate_ECrit`, psiUpstream, kxylemmax, c, d, pCrit)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_ECapacitance <- function(psi, psiPrev, PLCprev, V, fapo, c, d, pi0, eps, timestep) {
    .Call(`_medfate_ECapacitance`, psi, psiPrev, PLCprev, V, fapo, c, d, pi0, eps, timestep)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiVanGenuchten <- function(E, psiSoil, krhizomax, n, alpha, psiStep = -0.0001, psiMax = -10.0) {
    .Call(`_medfate_E2psiVanGenuchten`, E, psiSoil, krhizomax, n, alpha, psiStep, psiMax)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiTwoElements <- function(E, psiSoil, krhizomax, kxylemmax, n, alpha, c, d, psiCav = 0.0, psiStep = -0.0001, psiMax = -10.0) {
    .Call(`_medfate_E2psiTwoElements`, E, psiSoil, krhizomax, kxylemmax, n, alpha, c, d, psiCav, psiStep, psiMax)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiBelowground <- function(E, hydraulicNetwork, psiIni = as.numeric( c(0)), ntrial = 10L, psiTol = 0.0001, ETol = 0.0001) {
    .Call(`_medfate_E2psiBelowground`, E, hydraulicNetwork, psiIni, ntrial, psiTol, ETol)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiAboveground <- function(E, psiRootCrown, hydraulicNetwork) {
    .Call(`_medfate_E2psiAboveground`, E, psiRootCrown, hydraulicNetwork)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiFineRootLeaf <- function(E, psiFineRoot, hydraulicNetwork) {
    .Call(`_medfate_E2psiFineRootLeaf`, E, psiFineRoot, hydraulicNetwork)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiNetworkStem1 <- function(E, hydraulicNetwork, psiIni = as.numeric( c(0)), ntrial = 10L, psiTol = 0.0001, ETol = 0.0001) {
    .Call(`_medfate_E2psiNetworkStem1`, E, hydraulicNetwork, psiIni, ntrial, psiTol, ETol)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_E2psiNetwork <- function(E, hydraulicNetwork, psiIni = as.numeric( c(0)), ntrial = 10L, psiTol = 0.0001, ETol = 0.0001) {
    .Call(`_medfate_E2psiNetwork`, E, hydraulicNetwork, psiIni, ntrial, psiTol, ETol)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionOneXylem <- function(psiSoil, v, kstemmax, stemc, stemd, psiCav = 0.0, maxNsteps = 200L, dE = 0.01) {
    .Call(`_medfate_supplyFunctionOneXylem`, psiSoil, v, kstemmax, stemc, stemd, psiCav, maxNsteps, dE)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionTwoElements <- function(Emax, psiSoil, krhizomax, kxylemmax, n, alpha, c, d, psiCav = 0.0, dE = 0.1, psiMax = -10.0) {
    .Call(`_medfate_supplyFunctionTwoElements`, Emax, psiSoil, krhizomax, kxylemmax, n, alpha, c, d, psiCav, dE, psiMax)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionThreeElements <- function(Emax, psiSoil, krhizomax, kxylemmax, kleafmax, n, alpha, stemc, stemd, leafc, leafd, psiCav = 0.0, dE = 0.1, psiMax = -10.0) {
    .Call(`_medfate_supplyFunctionThreeElements`, Emax, psiSoil, krhizomax, kxylemmax, kleafmax, n, alpha, stemc, stemd, leafc, leafd, psiCav, dE, psiMax)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionBelowground <- function(hydraulicNetwork, minFlow = 0.0, maxNsteps = 400L, ntrial = 10L, psiTol = 0.0001, ETol = 0.0001, pCrit = 0.001) {
    .Call(`_medfate_supplyFunctionBelowground`, hydraulicNetwork, minFlow, maxNsteps, ntrial, psiTol, ETol, pCrit)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionAboveground <- function(Erootcrown, psiRootCrown, hydraulicNetwork) {
    .Call(`_medfate_supplyFunctionAboveground`, Erootcrown, psiRootCrown, hydraulicNetwork)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionFineRootLeaf <- function(psiFineRoot, hydraulicNetwork, minFlow = 0.0, maxNsteps = 400L, ETol = 0.0001, pCrit = 0.001) {
    .Call(`_medfate_supplyFunctionFineRootLeaf`, psiFineRoot, hydraulicNetwork, minFlow, maxNsteps, ETol, pCrit)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionNetworkStem1 <- function(hydraulicNetwork, minFlow = 0.0, maxNsteps = 400L, ntrial = 200L, psiTol = 0.0001, ETol = 0.0001, pCrit = 0.001) {
    .Call(`_medfate_supplyFunctionNetworkStem1`, hydraulicNetwork, minFlow, maxNsteps, ntrial, psiTol, ETol, pCrit)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_supplyFunctionNetwork <- function(hydraulicNetwork, minFlow = 0.0, maxNsteps = 400L, ntrial = 200L, psiTol = 0.0001, ETol = 0.0001, pCrit = 0.001) {
    .Call(`_medfate_supplyFunctionNetwork`, hydraulicNetwork, minFlow, maxNsteps, ntrial, psiTol, ETol, pCrit)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_regulatedPsiXylem <- function(E, psiUpstream, kxylemmax, c, d, psiStep = -0.01) {
    .Call(`_medfate_regulatedPsiXylem`, E, psiUpstream, kxylemmax, c, d, psiStep)
}

#' @rdname hydraulics_supplyfunctions
hydraulics_regulatedPsiTwoElements <- function(Emax, psiSoil, krhizomax, kxylemmax, n, alpha, c, d, dE = 0.1, psiMax = -10.0) {
    .Call(`_medfate_regulatedPsiTwoElements`, Emax, psiSoil, krhizomax, kxylemmax, n, alpha, c, d, dE, psiMax)
}

#' Scaling from conductivity to conductance
#' 
#' Functions used to scale from tissue conductivity to conductance of different elements of the continuum.
#' 
#' @param psiSoil Soil water potential (in MPa). A scalar or a vector depending on the function.
#' @param psiRhizo Water potential (in MPa) in the rhizosphere (root surface).
#' @param psiStem Water potential (in MPa) in the stem.
#' @param psiLeaf Water potential (in MPa) in the leaf.
#' @param PLCstem Percent loss of conductance (in \%) in the stem.
#' @param L Vector with the length of coarse roots (mm) for each soil layer.
#' @param V Vector with the proportion [0-1] of fine roots within each soil layer.
#' @param krhizomax Maximum rhizosphere hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param kleafmax Maximum leaf hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param kstemmax Maximum stem xylem hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param krootmax Maximum root xylem hydraulic conductance (defined as flow per leaf surface unit and per pressure drop).
#' @param psiStep Water potential precision (in MPa).
#' @param rootc,rootd Parameters of the Weibull function for roots (root xylem vulnerability curve).
#' @param stemc,stemd Parameters of the Weibull function for stems (stem xylem vulnerability curve).
#' @param leafc,leafd Parameters of the Weibull function for leaves (leaf vulnerability curve).
#' @param n,alpha Parameters of the Van Genuchten function (rhizosphere vulnerability curve).
#' @param averageResistancePercent Average (across water potential values) resistance percent of the rhizosphere, with respect to total resistance (rhizosphere + root xylem + stem xylem).
#' @param initialValue Initial value of rhizosphere conductance.
#' @param xylemConductivity Xylem conductivity as flow per length of conduit and pressure drop (in kg·m-1·s-1·MPa-1).
#' @param Al2As Leaf area to sapwood area (in m2·m-2).
#' @param height Plant height (in cm).
#' @param refheight Reference plant height of measurement of xylem conductivity (in cm).
#' @param taper A boolean flag to indicate correction by taper of xylem conduits (Christoffersen et al. 2017).
#' 
#' @details Details of the hydraulic model are given in the medfate book
#' 
#' @return
#' Values returned for each function are:
#' \itemize{
#'   \item{\code{hydraulics_maximumSoilPlantConductance}: The maximum soil-plant conductance, in the same units as the input segment conductances.}
#'   \item{\code{hydraulics_averageRhizosphereResistancePercent}: The average percentage of resistance due to the rhizosphere, calculated across water potential values.}
#'   \item{\code{hydraulics_findRhizosphereMaximumConductance}: The maximum rhizosphere conductance value given an average rhizosphere resistance and the vulnerability curves of rhizosphere, root and stem elements.}
#'   \item{\code{hydraulics_taperFactorSavage}: Taper factor according to Savage et al. (2010).}
#' }
#' 
#' @references
#' Christoffersen, B. O., M. Gloor, S. Fauset, N. M. Fyllas, D. R. Galbraith, T. R. Baker, L. Rowland, R. A. Fisher, O. J. Binks, S. A. Sevanto, C. Xu, S. Jansen, B. Choat, M. Mencuccini, N. G. McDowell, and P. Meir. 2016. Linking hydraulic traits to tropical forest function in a size-structured and trait-driven model (TFS v.1-Hydro). Geoscientific Model Development Discussions 9: 4227–4255.
#' 
#' Savage, V. M., L. P. Bentley, B. J. Enquist, J. S. Sperry, D. D. Smith, P. B. Reich, and E. I. von Allmen. 2010. Hydraulic trade-offs and space filling enable better predictions of vascular structure and function in plants. Proceedings of the National Academy of Sciences of the United States of America 107:22722–7.
#' 
#' Olson, M.E., Anfodillo, T., Rosell, J.A., Petit, G., Crivellaro, A., Isnard, S., \enc{León-Gómez}{Leon-Gomez}, C., \enc{Alvarado-Cárdenas}{Alvarado-Cardenas}, L.O., and Castorena, M. 2014. Universal hydraulics of the flowering plants: Vessel diameter scales with stem length across angiosperm lineages, habits and climates. Ecology Letters 17: 988–997.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{hydraulics_psi2K}}, \code{\link{hydraulics_supplyFunctionPlot}}, \code{\link{spwb}}, \code{\link{soil}}
#' 
#' @name hydraulics_scalingconductance
hydraulics_maximumSoilPlantConductance <- function(krhizomax, krootmax, kstemmax, kleafmax) {
    .Call(`_medfate_maximumSoilPlantConductance`, krhizomax, krootmax, kstemmax, kleafmax)
}

#' @rdname hydraulics_scalingconductance
hydraulics_soilPlantResistances <- function(psiSoil, psiRhizo, psiStem, PLCstem, psiLeaf, krhizomax, n, alpha, krootmax, rootc, rootd, kstemmax, stemc, stemd, kleafmax, leafc, leafd) {
    .Call(`_medfate_soilPlantResistances`, psiSoil, psiRhizo, psiStem, PLCstem, psiLeaf, krhizomax, n, alpha, krootmax, rootc, rootd, kstemmax, stemc, stemd, kleafmax, leafc, leafd)
}

#' @rdname hydraulics_scalingconductance
hydraulics_averageRhizosphereResistancePercent <- function(krhizomax, n, alpha, krootmax, rootc, rootd, kstemmax, stemc, stemd, kleafmax, leafc, leafd, psiStep = -0.01) {
    .Call(`_medfate_averageRhizosphereResistancePercent`, krhizomax, n, alpha, krootmax, rootc, rootd, kstemmax, stemc, stemd, kleafmax, leafc, leafd, psiStep)
}

#' @rdname hydraulics_scalingconductance
hydraulics_findRhizosphereMaximumConductance <- function(averageResistancePercent, n, alpha, krootmax, rootc, rootd, kstemmax, stemc, stemd, kleafmax, leafc, leafd, initialValue = 0.0) {
    .Call(`_medfate_findRhizosphereMaximumConductance`, averageResistancePercent, n, alpha, krootmax, rootc, rootd, kstemmax, stemc, stemd, kleafmax, leafc, leafd, initialValue)
}

#' @rdname hydraulics_scalingconductance
hydraulics_taperFactorSavage <- function(height) {
    .Call(`_medfate_taperFactorSavage`, height)
}

#' @rdname hydraulics_scalingconductance
hydraulics_terminalConduitRadius <- function(height) {
    .Call(`_medfate_terminalConduitRadius`, height)
}

#' @rdname hydraulics_scalingconductance
hydraulics_referenceConductivityHeightFactor <- function(refheight, height) {
    .Call(`_medfate_referenceConductivityHeightFactor`, refheight, height)
}

#' @rdname hydraulics_scalingconductance
hydraulics_maximumStemHydraulicConductance <- function(xylemConductivity, refheight, Al2As, height, taper = FALSE) {
    .Call(`_medfate_maximumStemHydraulicConductance`, xylemConductivity, refheight, Al2As, height, taper)
}

#' @rdname hydraulics_scalingconductance
hydraulics_rootxylemConductanceProportions <- function(L, V) {
    .Call(`_medfate_rootxylemConductanceProportions`, L, V)
}

#' @param doy Day of the year.
#' @param pet Potential evapotranspiration for a given day (mm).
#' @param prec Precipitation for a given day (mm).
#' @param Rconv,Rsyn Rainfall rate for convective storms and synoptic storms, respectively, in mm/h.
#' 
#' @rdname hydrology_interception
hydrology_erFactor <- function(doy, pet, prec, Rconv = 5.6, Rsyn = 1.5) {
    .Call(`_medfate_erFactor`, doy, pet, prec, Rconv, Rsyn)
}

.hydrology_interceptionGashDay <- function(Precipitation, Cm, p, ER = 0.05) {
    .Call(`_medfate_interceptionGashDay`, Precipitation, Cm, p, ER)
}

#' @rdname hydrology_soil
#' 
#' @param DEF Water deficit in the (topsoil) layer.
#' @param PETs Potential evapotranspiration at the soil surface.
#' @param Gsoil Gamma parameter (maximum daily evaporation).
#' 
hydrology_soilEvaporationAmount <- function(DEF, PETs, Gsoil) {
    .Call(`_medfate_soilEvaporationAmount`, DEF, PETs, Gsoil)
}

#' @rdname hydrology_soil
#' 
#' @param soil An object of class \code{\link{soil}}.
#' @param soilFunctions Soil water retention curve and conductivity functions, either 'SX' (for Saxton) or 'VG' (for Van Genuchten).
#' @param pet Potential evapotranspiration for a given day (mm)
#' @param LgroundSWR Percentage of short-wave radiation (SWR) reaching the ground.
#' @param modifySoil Boolean flag to indicate that the input \code{soil} object should be modified during the simulation.
#' 
hydrology_soilEvaporation <- function(soil, soilFunctions, pet, LgroundSWR, modifySoil = TRUE) {
    .Call(`_medfate_soilEvaporation`, soil, soilFunctions, pet, LgroundSWR, modifySoil)
}

#' @rdname hydrology_soil
#' @param LherbSWR Percentage of short-wave radiation (SWR) reaching the herbaceous layer.
#' @param herbLAI Leaf area index of the herbaceous layer.
hydrology_herbaceousTranspiration <- function(pet, LherbSWR, herbLAI, soil, soilFunctions, modifySoil = TRUE) {
    .Call(`_medfate_herbaceousTranspiration`, pet, LherbSWR, herbLAI, soil, soilFunctions, modifySoil)
}

.hydrology_infiltrationAmount <- function(input, Ssoil) {
    .Call(`_medfate_infiltrationAmount`, input, Ssoil)
}

#' @rdname hydrology_soil
#' 
#' @param I Soil infiltration (in mm of water).
#' @param dVec Width of soil layers (in mm).
#' @param macro Macroporosity of soil layers (in \%).
#' @param a,b Parameters of the extinction function used for water infitration.
#' 
hydrology_infiltrationRepartition <- function(I, dVec, macro, a = -0.005, b = 3.0) {
    .Call(`_medfate_infiltrationRepartition`, I, dVec, macro, a, b)
}

#' @rdname hydrology_soil
#' 
#' @param tday Average day temperature (ºC).
#' @param rad Solar radiation (in MJ/m2/day).
#' @param elevation Altitude above sea level (m).
#' 
hydrology_snowMelt <- function(tday, rad, LgroundSWR, elevation) {
    .Call(`_medfate_snowMelt`, tday, rad, LgroundSWR, elevation)
}

#' Soil vertical inputs
#' 
#' High-level functions for hydrological processes. Function \code{hydrology_soilWaterInputs} performs 
#' canopy water interception and snow accumulation/melt. Function \code{hydrology_soilInfiltrationPercolation} 
#' performs soil infiltration and percolation from the input given by the previous function.
#' 
#' @param soil A list containing the description of the soil (see \code{\link{soil}}).
#' @param soilFunctions Soil water retention curve and conductivity functions, either 'SX' (for Saxton) or 'VG' (for Van Genuchten).
#' @param prec Precipitation for a given day (mm)
#' @param er The ratio of evaporation rate to rainfall rate.
#' @param tday Average day temperature (ºC).
#' @param rad Solar radiation (in MJ/m2/day).
#' @param elevation Altitude above sea level (m).
#' @param Cm Canopy water storage capacity.
#' @param LgroundPAR Percentage of photosynthetically-acvive radiation (PAR) reaching the ground.
#' @param LgroundSWR Percentage of short-wave radiation (SWR) reaching the ground.
#' @param runon Surface water amount running on the target area from upslope (in mm).
#' @param snowpack Boolean flag to indicate the simulation of snow accumulation and melting.
#' @param modifySoil Boolean flag to indicate that the input \code{soil} object should be modified during the simulation.
#' 
#' @details 
#' The function simulates different vertical hydrological processes, which are described separately in other functions. 
#' If \code{modifySoil = TRUE} the function will modify the \code{soil} object (including both soil moisture and 
#' the snowpack on its surface) as a result of simulating hydrological processes.
#' 
#' @return 
#' Function \code{hydrology_soilWaterInputs} returns a named vector with the following elements, all in mm:
#' \item{Rain}{Precipitation as rainfall.}
#' \item{Snow}{Precipitation as snow.}
#' \item{Interception}{Rainfall water intercepted by the canopy and evaporated.}
#' \item{NetRain}{Rainfall reaching the ground.}
#' \item{Snowmelt}{Snow melted during the day, and added to the water infiltrated.}
#' \item{Runon}{Surface water amount running on the target area from upslope.}
#' \item{Input}{Total soil input, including runon, snowmelt and net rain.}
#' 
#' Function \code{hydrology_soilInfiltrationPercolation} returns a named vector with the following elements, all in mm:
#' \item{Infiltration}{Water infiltrated into the soil (i.e. throughfall + runon + snowmelt - runoff).}
#' \item{Runoff}{Surface water leaving the target area.}
#' \item{DeepDrainage}{Water leaving the target soil towards the water table.}
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{spwb_day}}, \code{\link{hydrology_rainInterception}}, \code{\link{hydrology_soilEvaporation}}
#' 
#' @name hydrology_verticalInputs
hydrology_soilWaterInputs <- function(soil, soilFunctions, prec, er, tday, rad, elevation, Cm, LgroundPAR, LgroundSWR, runon = 0.0, snowpack = TRUE, modifySoil = TRUE) {
    .Call(`_medfate_soilWaterInputs`, soil, soilFunctions, prec, er, tday, rad, elevation, Cm, LgroundPAR, LgroundSWR, runon, snowpack, modifySoil)
}

#' @rdname hydrology_verticalInputs
#' 
#' @param waterInput Soil water input for a given day (mm).
#' @param rockyLayerDrainage Boolean flag to indicate the simulation of drainage from rocky layers (> 95\% of rocks).
#' 
hydrology_soilInfiltrationPercolation <- function(soil, soilFunctions, waterInput, rockyLayerDrainage = TRUE, modifySoil = TRUE) {
    .Call(`_medfate_soilInfiltrationPercolation`, soil, soilFunctions, waterInput, rockyLayerDrainage, modifySoil)
}

.incgam <- function(a, x) {
    .Call(`_medfate_incgam`, a, x)
}

.invincgam <- function(a, p, q) {
    .Call(`_medfate_invincgam`, a, p, q)
}

.gammds <- function(x, p) {
    .Call(`_medfate_gammds`, x, p)
}

#' Sureau-ECOS inner functions for testing only
#' 
#' Function \code{initCochardNetworks} initializes hydraulic networks for all plant cohorts in x
#' Function \code{semi_implicit_integration} updates water potentials and cavitation across the hydraulic network
#' 
#' @param x An object of class \code{\link{spwbInput}} or \code{\link{growthInput}} created using \code{transpirationMode = "Cochard"}.
#'  
#' @return Function \code{initCochardNetworks} returns a vector of length equal to the number of cohorts. Each element is a list with Sureau-ECOS parameters.
#' Function \code{semi_implicit_integration} does not return anything, but modifies input parameter \code{network}.
#' 
#' @author
#' \itemize{
#'   \item{Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF}
#'   \item{Nicolas Martin-StPaul, URFM-INRAE}
#' }
#' 
#' @references
#' Ruffault J, Pimont F, Cochard H, Dupuy JL, Martin-StPaul N (2022) 
#' SurEau-Ecos v2.0: a trait-based plant hydraulics model for simulations of plant water status and drought-induced mortality at the ecosystem level.
#' Geoscientific Model Development 15, 5593-5626 (doi:10.5194/gmd-15-5593-2022).
#' 
#' 
#' @seealso  \code{\link{spwb}}
#' 
#' @name sureau_ecos
initCochardNetworks <- function(x) {
    .Call(`_medfate_initCochardNetworks`, x)
}

#' @rdname sureau_ecos
#' @param network A hydraulic network element of the list returned by \code{initCochardNetworks}
#' @param dt Smallest time step (seconds)
#' @param opt Option flag vector
semi_implicit_integration <- function(network, dt, opt) {
    invisible(.Call(`_medfate_semi_implicit_integration`, network, dt, opt))
}

#' Stomatal regulation
#' 
#' Set of high-level functions used in the calculation of stomatal conductance and transpiration. 
#' Function \code{transp_profitMaximization} calculates gain and cost functions, 
#' as well as profit maximization from supply and photosynthesis input functions. 
#' Function \code{transp_stomatalRegulationPlot} produces a plot with the cohort supply functions against water potential 
#' and a plot with the cohort photosynthesis functions against water potential, both with the maximum profit values indicated.
#' 
#' @param supplyFunction Water supply function (see \code{\link{hydraulics_supplyFunctionNetwork}}).
#' @param photosynthesisFunction Function returned by \code{photo_photosynthesisFunction()}.
#' @param Gswmin,Gswmax Minimum and maximum stomatal conductance to water vapour (mol·m-2·s-1).
#' 
#' @return
#' Function \code{transp_profitMaximization} returns a list with the following elements:
#' \itemize{
#'   \item{\code{Cost}: Cost function [0-1].}
#'   \item{\code{Gain}: Gain function [0-1].}
#'   \item{\code{Profit}: Profit function [0-1].}
#'   \item{\code{iMaxProfit}: Index corresponding to maximum profit (starting from 0).}
#' }
#' 
#' @references
#' Sperry, J. S., M. D. Venturas, W. R. L. Anderegg, M. Mencuccini, D. S. Mackay, Y. Wang, and D. M. Love. 2017. Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost. Plant Cell and Environment 40, 816-830 (doi: 10.1111/pce.12852).
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{transp_transpirationSperry}}, \code{\link{hydraulics_supplyFunctionNetwork}}, \code{\link{biophysics_leafTemperature}}, \code{\link{photo_photosynthesis}}, \code{\link{spwb_day}}, \code{\link{plot.spwb_day}}
#' 
#' @examples
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Initialize soil with default soil params (4 layers)
#' examplesoil = soil(defaultSoilParams(4))
#' 
#' #Initialize control parameters
#' control = defaultControl(transpirationMode="Sperry")
#' 
#' #Initialize input
#' x2 = forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' # Stomatal VPD curve and chosen value for the 12th time step at day 100
#' transp_stomatalRegulationPlot(x2, examplemeteo, day=100, timestep = 12,
#'                               latitude = 41.82592, elevation = 100, type="VPD")
#'  
#' @name transp_stomatalregulation
transp_profitMaximization <- function(supplyFunction, photosynthesisFunction, Gswmin, Gswmax) {
    .Call(`_medfate_profitMaximization`, supplyFunction, photosynthesisFunction, Gswmin, Gswmax)
}

.parcohort <- function(SP, H, CR, LAI, SpParams) {
    .Call(`_medfate_parcohort`, SP, H, CR, LAI, SpParams)
}

#' Radiation transfer functions
#' 
#' Functions \code{light_layerIrradianceFraction} and \code{light_layerIrradianceFractionBottomUp} calculate 
#' the fraction of above-canopy irradiance (and the soil irradiance, respectively) reaching each vegetation layer. 
#' Function \code{light_layerSunlitFraction} calculates the proportion of sunlit leaves in each vegetation layer. 
#' Function \code{light_cohortSunlitShadeAbsorbedRadiation} calculates the amount of radiation absorved 
#' by cohort and vegetation layers, while differentiating between sunlit and shade leaves.
#' 
#' @param LAIme A numeric matrix of live expanded LAI values per vegetation layer (row) and cohort (column).
#' @param LAImd A numeric matrix of dead LAI values per vegetation layer (row) and cohort (column).
#' @param LAImx A numeric matrix of maximum LAI values per vegetation layer (row) and cohort (column).
#' @param k A vector of light extinction coefficients.
#' @param kb A vector of direct light extinction coefficients.
#' @param kd A vector of diffuse light extinction coefficients.
#' @param Ib0 Above-canopy direct incident radiation.
#' @param Id0 Above-canopy diffuse incident radiation.
#' @param Ibf Fraction of above-canopy direct radiation reaching each vegetation layer.
#' @param Idf Fraction of above-canopy diffuse radiation reaching each vegetation layer.
#' @param alpha A vecfor of leaf absorbance by species.
#' @param beta Solar elevation (in radians).
#' @param gamma Vector of canopy reflectance values.
#' @param kPAR A vector of visible light extinction coefficients for each cohort.
#' @param alphaSWR A vecfor of hort-wave absorbance coefficients for each cohort.
#' @param gammaSWR A vector of short-wave reflectance coefficients (albedo) for each cohort.
#' @param ddd A dataframe with direct and diffuse radiation for different subdaily time steps (see function \code{radiation_directDiffuseDay} in package meteoland).
#' @param ntimesteps Number of subdaily time steps.
#' @param trunkExtinctionFraction Fraction of extinction due to trunks (for winter deciduous forests).
#' @param LWRatm Atmospheric downward long-wave radiation (W/m2).
#' @param Tsoil Soil temperature (Celsius).
#' @param Tair Canopy layer air temperature vector (Celsius).
#' @param x An object of class \code{\link{forest}}
#' @param SpParams A data frame with species parameters (see \code{\link{SpParamsMED}}).
#' @param z A numeric vector with height values.
#' @param gdd Growth degree days.
#' 
#' @details
#' Functions for short-wave radiation are adapted from Anten & Bastiaans (2016), 
#' whereas long-wave radiation balance follows Flerchinger et al. (2009). 
#' Vegetation layers are assumed to be ordered from bottom to top.
#' 
#' @return
#' Functions \code{light_layerIrradianceFraction}, \code{light_layerIrradianceFractionBottomUp}  and \code{light_layerSunlitFraction} 
#' return a numeric vector of length equal to the number of vegetation layers. 
#' 
#' Function \code{light_cohortSunlitShadeAbsorbedRadiation} returns a list with 
#' two elements (matrices): \code{I_sunlit} and \code{I_shade}.
#' 
#' @references
#' Anten, N.P.R., Bastiaans, L., 2016. The use of canopy models to analyze light competition among plants, in: Hikosaka, K., Niinemets, U., Anten, N.P.R. (Eds.), Canopy Photosynthesis: From Basics to Application. Springer, pp. 379–398.
#' 
#' Flerchinger, G. N., Xiao, W., Sauer, T. J., Yu, Q. 2009. Simulation of within-canopy radiation exchange. NJAS - Wageningen Journal of Life Sciences 57 (1): 5–15. https://doi.org/10.1016/j.njas.2009.07.004.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso  \code{\link{spwb}}
#' 
#' @examples
#' LAI <- 2
#' nlayer <- 10
#' LAIlayerlive <- matrix(rep(LAI/nlayer,nlayer),nlayer,1)
#' LAIlayerdead <- matrix(0,nlayer,1)
#' kb <- 0.8
#' kd_PAR <- 0.5
#' kd_SWR <- kd_PAR/1.35
#' alpha_PAR <- 0.9
#' gamma_PAR <- 0.04
#' gamma_SWR <- 0.05
#' alpha_SWR <- 0.7
#' 
#' Ibfpar <- light_layerIrradianceFraction(LAIlayerlive,LAIlayerdead,LAIlayerlive,kb, alpha_PAR)
#' Idfpar <- light_layerIrradianceFraction(LAIlayerlive,LAIlayerdead,LAIlayerlive,kd_PAR, alpha_PAR)
#' Ibfswr <- light_layerIrradianceFraction(LAIlayerlive,LAIlayerdead,LAIlayerlive,kb, alpha_SWR)
#' Idfswr <- light_layerIrradianceFraction(LAIlayerlive,LAIlayerdead,LAIlayerlive,kd_SWR, alpha_SWR)
#' fsunlit <- light_layerSunlitFraction(LAIlayerlive, LAIlayerdead, kb)
#' SHarea <- (1-fsunlit)*LAIlayerlive[,1] 
#' SLarea <- fsunlit*LAIlayerlive[,1] 
#' 
#' oldpar <- par(mar=c(4,4,1,1), mfrow=c(1,2))
#' plot(Ibfpar*100, 1:nlayer,type="l", ylab="Layer", 
#'      xlab="Percentage of irradiance", xlim=c(0,100), ylim=c(1,nlayer), col="dark green")
#' lines(Idfpar*100, 1:nlayer, col="dark green", lty=2)
#' lines(Ibfswr*100, 1:nlayer, col="red")
#' lines(Idfswr*100, 1:nlayer, col="red", lty=2)
#'   
#' plot(fsunlit*100, 1:nlayer,type="l", ylab="Layer", 
#'      xlab="Percentage of leaves", xlim=c(0,100), ylim=c(1,nlayer))
#' lines((1-fsunlit)*100, 1:nlayer, lty=2)
#' par(oldpar)  
#'   
#' solarElevation <- 0.67
#' SWR_direct <- 1100
#' SWR_diffuse <- 300
#' PAR_direct <- 550
#' PAR_diffuse <- 150
#' 
#' abs_PAR <- light_cohortSunlitShadeAbsorbedRadiation(PAR_direct, PAR_diffuse,
#'                         Ibfpar, Idfpar, beta = solarElevation, 
#'                         LAIlayerlive, LAIlayerdead, kb, kd_PAR, alpha_PAR, gamma_PAR)
#' abs_SWR <- light_cohortSunlitShadeAbsorbedRadiation(SWR_direct, SWR_diffuse,
#'                          Ibfswr, Idfswr, beta = solarElevation, 
#'                          LAIlayerlive, LAIlayerdead, kb, kd_SWR, alpha_SWR, gamma_SWR)
#' oldpar <- par(mar=c(4,4,1,1), mfrow=c(1,2))
#' absRadSL <- abs_SWR$I_sunlit[,1]
#' absRadSH <- abs_SWR$I_shade[,1]
#' lambda <- 546.6507
#' QSL <- abs_PAR$I_sunlit[,1]*lambda*0.836*0.01
#' QSH <- abs_PAR$I_shade[,1]*lambda*0.836*0.01
#' plot(QSL, 1:nlayer,type="l", ylab="Layer", 
#'    xlab="Absorbed PAR quantum flux per leaf area", ylim=c(1,nlayer), col="dark green", 
#'    xlim=c(0,max(QSL)))
#' lines(QSH, 1:nlayer, col="dark green", lty=2)
#' plot(absRadSL, 1:nlayer,type="l", ylab="Layer", 
#'    xlab="Absorbed SWR per leaf area (W/m2)", ylim=c(1,nlayer), col="red", 
#'    xlim=c(0,max(absRadSL)))
#' lines(absRadSH, 1:nlayer, col="red", lty=2)
#' par(oldpar)
#'   
#' @name light
light_PARcohort <- function(x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_PARcohort`, x, SpParams, gdd)
}

.parheight <- function(z, x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_parheight`, z, x, SpParams, gdd)
}

#' @rdname light
light_PARground <- function(x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_PARground`, x, SpParams, gdd)
}

.swrheight <- function(z, x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_swrheight`, z, x, SpParams, gdd)
}

#' @rdname light
light_SWRground <- function(x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_SWRground`, x, SpParams, gdd)
}

.parExtinctionProfile <- function(z, x, SpParams, gdd = NA_real_, includeHerbs = FALSE) {
    .Call(`_medfate_parExtinctionProfile`, z, x, SpParams, gdd, includeHerbs)
}

.swrExtinctionProfile <- function(z, x, SpParams, gdd = NA_real_, includeHerbs = FALSE) {
    .Call(`_medfate_swrExtinctionProfile`, z, x, SpParams, gdd, includeHerbs)
}

#' @rdname light
light_cohortAbsorbedSWRFraction <- function(z, x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_cohortAbsorbedSWRFraction`, z, x, SpParams, gdd)
}

#' @rdname light
light_layerIrradianceFraction <- function(LAIme, LAImd, LAImx, k, alpha, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_layerIrradianceFraction`, LAIme, LAImd, LAImx, k, alpha, trunkExtinctionFraction)
}

#' @rdname light
light_layerIrradianceFractionBottomUp <- function(LAIme, LAImd, LAImx, k, alpha, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_layerIrradianceFractionBottomUp`, LAIme, LAImd, LAImx, k, alpha, trunkExtinctionFraction)
}

#' @rdname light
light_cohortSunlitShadeAbsorbedRadiation <- function(Ib0, Id0, Ibf, Idf, beta, LAIme, LAImd, kb, kd, alpha, gamma) {
    .Call(`_medfate_cohortSunlitShadeAbsorbedRadiation`, Ib0, Id0, Ibf, Idf, beta, LAIme, LAImd, kb, kd, alpha, gamma)
}

#' @rdname light
light_layerSunlitFraction <- function(LAIme, LAImd, kb) {
    .Call(`_medfate_layerSunlitFraction`, LAIme, LAImd, kb)
}

#' @rdname light
light_instantaneousLightExtinctionAbsortion <- function(LAIme, LAImd, LAImx, kPAR, alphaSWR, gammaSWR, ddd, ntimesteps = 24L, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_instantaneousLightExtinctionAbsortion`, LAIme, LAImd, LAImx, kPAR, alphaSWR, gammaSWR, ddd, ntimesteps, trunkExtinctionFraction)
}

#' @rdname light
light_longwaveRadiationSHAW <- function(LAIme, LAImd, LAImx, LWRatm, Tsoil, Tair, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_longwaveRadiationSHAW`, LAIme, LAImd, LAImx, LWRatm, Tsoil, Tair, trunkExtinctionFraction)
}

.paramsBelow <- function(above, Z50, Z95, soil, paramsAnatomydf, paramsTranspirationdf, control) {
    .Call(`_medfate_paramsBelow`, above, Z50, Z95, soil, paramsAnatomydf, paramsTranspirationdf, control)
}

.spwbInput <- function(above, Z50, Z95, soil, FCCSprops, SpParams, control) {
    .Call(`_medfate_spwbInput`, above, Z50, Z95, soil, FCCSprops, SpParams, control)
}

.growthInput <- function(above, Z50, Z95, soil, FCCSprops, SpParams, control) {
    .Call(`_medfate_growthInput`, above, Z50, Z95, soil, FCCSprops, SpParams, control)
}

.cloneInput <- function(input) {
    .Call(`_medfate_cloneInput`, input)
}

#' Input for simulation models
#'
#' Functions \code{forest2spwbInput} and \code{forest2growthInput} take an object of class \code{\link{forest}} 
#' and create input objects for simulation functions \code{\link{spwb}} (or \code{\link{pwb}}) and \code{\link{growth}}, respectively. 
#' Function \code{forest2aboveground} calculates aboveground variables such as leaf area index. 
#' Function \code{forest2belowground} calculates belowground variables such as fine root distribution.
#' 
#' @param x An object of class \code{\link{forest}}.
#' @param SpParams A data frame with species parameters (see \code{\link{SpParamsDefinition}} and \code{\link{SpParamsMED}}).
#' @param gdd Growth degree days to account for leaf phenology effects (in Celsius). This should be left \code{NA} in most applications.
#' @param loading A logical flag to indicate that fuel loading should be included (for fire hazard calculations). 
#' @param soil An object of class \code{\link{soil}}.
#' @param control A list with default control parameters (see \code{\link{defaultControl}}).
#' 
#' @details
#' Function \code{forest2aboveground} extract height and species identity from plant cohorts of \code{x}, 
#' and calculate leaf area index and crown ratio. Functions \code{forest2spwbInput} and \code{forest2growthInput} also calculate the distribution of fine roots 
#' across soil, and finds parameter values for each plant cohort according to the parameters of its species as specified in \code{SpParams}. 
#' If \code{control$transpirationMode = "Sperry"} or \code{control$transpirationMode = "Cochard"},
#' the \code{forest2spwbInput} and \code{forest2growthInput} also estimate the maximum conductance of rhizosphere, root xylem and stem xylem elements.
#' 
#' @return 
#' Function \code{forest2aboveground()} returns a data frame with the following columns (rows are identified as specified by function \code{\link{plant_ID}}):
#' \itemize{
#'   \item{\code{SP}: Species identity (an integer) (first species is 0).}
#'   \item{\code{N}: Cohort density (ind/ha) (see function \code{\link{plant_density}}).}
#'   \item{\code{DBH}: Tree diameter at breast height (cm).}
#'   \item{\code{H}: Plant total height (cm).}
#'   \item{\code{CR}: Crown ratio (crown length to total height) (between 0 and 1).}
#'   \item{\code{LAI_live}: Live leaf area index (m2/m2) (one-side leaf area relative to plot area), includes leaves in winter dormant buds.}
#'   \item{\code{LAI_expanded}: Leaf area index of expanded leaves (m2/m2) (one-side leaf area relative to plot area).}
#'   \item{\code{LAI_dead}: Dead leaf area index (m2/m2) (one-side leaf area relative to plot area).}
#'   \item{\code{Loading}: Fine fuel loading (kg/m2), only if \code{loading = TRUE}.}
#' }
#' 
#' Function \code{forest2spwbInput()} returns a list of class \code{spwbInput} with the following elements (rows of data frames are identified as specified by function \code{\link{plant_ID}}):
#'   \itemize{
#'     \item{\code{control}: List with control parameters (see \code{\link{defaultControl}}).}
#'     \item{\code{canopy}: A list of stand-level state variables.}
#'     \item{\code{cohorts}: A data frame with cohort information, with columns \code{SP} and \code{Name}.}
#'     \item{\code{above}: A data frame with columns  \code{H}, \code{CR} and \code{LAI} (see function \code{forest2aboveground}).}
#'     \item{\code{below}: A data frame with columns \code{Z50}, \code{Z95}.  If \code{control$transpirationMode = "Sperry"} additional columns are \code{fineRootBiomass} and \code{coarseRootSoilVolume}.}
#'     \item{\code{belowLayers}: A list. If \code{control$transpirationMode = "Granier"} it contains elements: 
#'       \itemize{
#'         \item{\code{V}: A matrix with the proportion of fine roots of each cohort (in rows) in each soil layer (in columns).}
#'         \item{\code{L}: A matrix with the length of coarse roots of each cohort (in rows) in each soil layer (in columns).}
#'         \item{\code{Wpool}: A matrix with the soil moisture relative to field capacity around the rhizosphere of each cohort (in rows) in each soil layer (in columns).}
#'       }
#'       If \code{control$transpirationMode = "Sperry"} or \code{control$transpirationMode = "Cochard"} there are the following additional elements:
#'       \itemize{
#'         \item{\code{VGrhizo_kmax}: A matrix with maximum rhizosphere conductance values of each cohort (in rows) in each soil layer (in columns).}
#'         \item{\code{VGroot_kmax}: A matrix with maximum root xylem conductance values of each cohort (in rows) in each soil layer (in columns).}
#'         \item{\code{RhizoPsi}: A matrix with the water potential around the rhizosphere of each cohort (in rows) in each soil layer (in columns).}
#'       }
#'     }
#'     \item{\code{paramsPhenology}: A data frame with leaf phenology parameters:
#'       \itemize{
#'         \item{\code{PhenologyType}: Leaf phenology type.}
#'         \item{\code{LeafDuration}: Leaf duration (in years).}
#'         \item{\code{Sgdd}: Degree days needed for leaf budburst (for winter decideous species).}
#'         \item{\code{Tbgdd}: Base temperature for the calculation of degree days to leaf budburst.}
#'         \item{\code{Ssen}: Degree days corresponding to leaf senescence.}
#'         \item{\code{Phsen}: Photoperiod corresponding to start counting senescence degree-days.}
#'         \item{\code{Tbsen}: Base temperature for the calculation of degree days to leaf senescence.}
#'       }
#'     }
#'     \item{\code{paramsAnatomy}: A data frame with plant anatomy parameters for each cohort:
#'       \itemize{
#'         \item{\code{Hmax}: Maximum plant height (cm).}
#'         \item{\code{Hmed}: Median plant height (cm).}
#'         \item{\code{Al2As}: Leaf area to sapwood area ratio (in m2·m-2).}
#'         \item{\code{Ar2Al}: Fine root area to leaf area ratio (in m2·m-2).}
#'         \item{\code{SLA}: Specific leaf area (mm2/mg = m2/kg).}
#'         \item{\code{LeafWidth}: Leaf width (in cm).}
#'         \item{\code{LeafDensity}: Density of leaf tissue (dry weight over volume).}
#'         \item{\code{WoodDensity}: Density of wood tissue (dry weight over volume).}
#'         \item{\code{FineRootDensity}: Density of fine root tissue (dry weight over volume).}
#'         \item{\code{SRL}: Specific Root length (cm·g-1).}
#'         \item{\code{RLD}: Root length density (cm·cm-3).}
#'         \item{\code{r635}: Ratio between the weight of leaves plus branches and the weight of leaves alone for branches of 6.35 mm.}
#'       }
#'     }
#'     \item{\code{paramsInterception}: A data frame with rain interception and light extinction parameters for each cohort:
#'       \itemize{
#'         \item{\code{kPAR}: PAR extinction coefficient.}
#'         \item{\code{g}: Canopy water retention capacity per LAI unit (mm/LAI).}
#'       }
#'     If \code{control$transpirationMode = "Sperry"} or \code{control$transpirationMode = "Cochard"} additional columns are:
#'       \itemize{
#'         \item{\code{gammaSWR}: Reflectance (albedo) coefficient for SWR .}
#'         \item{\code{alphaSWR}: Absorbance coefficient for SWR .}
#'       }
#'     }
#'     \item{\code{paramsTranspiration}: A data frame with parameters for transpiration and photosynthesis. If \code{control$transpirationMode = "Granier"}, columns are:
#'       \itemize{
#'         \item{\code{Gswmin}: Minimum stomatal conductance to water vapor (in mol H2O·m-2·s-1).}
#'         \item{\code{Tmax_LAI}: Coefficient relating LAI with the ratio of maximum transpiration over potential evapotranspiration.}
#'         \item{\code{Tmax_LAIsq}: Coefficient relating squared LAI with the ratio of maximum transpiration over potential evapotranspiration.}
#'         \item{\code{Psi_Extract}: Water potential corresponding to 50\% relative transpiration (in MPa).}
#'         \item{\code{Exp_Extract}: Parameter of the Weibull function regulating transpiration reduction.}
#'         \item{\code{VCstem_c}, \code{VCstem_d}: Parameters of the stem xylem vulnerability curve.}
#'         \item{\code{WUE}: Daily water use efficiency (gross photosynthesis over transpiration) under no light, water or CO2 limitations and VPD = 1kPa (g C/mm water).}
#'         \item{\code{WUE_par}: Coefficient regulating the influence of \% PAR on gross photosynthesis.}
#'         \item{\code{WUE_par}: Coefficient regulating the influence of atmospheric CO2 concentration on gross photosynthesis.}
#'         \item{\code{WUE_par}: Coefficient regulating the influence of vapor pressure deficit (VPD) on gross photosynthesis.}
#'       }
#'      If \code{control$transpirationMode = "Sperry"} columns are:
#'       \itemize{
#'         \item{\code{Gswmin}: Minimum stomatal conductance to water vapor (in mol H2O·m-2·s-1).}
#'         \item{\code{Gswmax}: Maximum stomatal conductance to water vapor (in mol H2O·m-2·s-1).}
#'         \item{\code{Vmax298}: Maximum Rubisco carboxilation rate at 25ºC (in micromol CO2·s-1·m-2).}
#'         \item{\code{Jmax298}: Maximum rate of electron transport at 25ºC (in micromol photons·s-1·m-2).}
#'         \item{\code{Kmax_stemxylem}: Sapwood-specific hydraulic conductivity of stem xylem (in kg H2O·s-1·m-2).}
#'         \item{\code{Kmax_rootxylem}: Sapwood-specific hydraulic conductivity of root xylem (in kg H2O·s-1·m-2).}
#'         \item{\code{VCleaf_kmax}: Maximum leaf hydraulic conductance.}
#'         \item{\code{VCleaf_c}, \code{VCleaf_d}: Parameters of the leaf vulnerability curve.}
#'         \item{\code{VCstem_kmax}: Maximum stem xylem conductance.}
#'         \item{\code{VCstem_c}, \code{VCstem_d}: Parameters of the stem xylem vulnerability curve.}
#'         \item{\code{VCroot_c}, \code{VCroot_d}: Parameters of the root xylem vulnerability curve.}
#'         \item{\code{Plant_kmax}: Maximum whole-plant conductance.}
#'       }
#'       If \code{control$transpirationMode = "Cochard"} columns are:
#'       \itemize{
#'         \item{\code{Gswmin}: Minimum stomatal conductance to water vapor (in mol H2O·m-2·s-1).}
#'         \item{\code{Gswmax}: Maximum stomatal conductance to water vapor (in mol H2O·m-2·s-1).}
#'         \item{\code{Vmax298}: Maximum Rubisco carboxilation rate at 25ºC (in micromol CO2·s-1·m-2).}
#'         \item{\code{Jmax298}: Maximum rate of electron transport at 25ºC (in micromol photons·s-1·m-2).}
#'         \item{\code{Kmax_stemxylem}: Sapwood-specific hydraulic conductivity of stem xylem (in kg H2O·s-1·m-2).}
#'         \item{\code{Kmax_rootxylem}: Sapwood-specific hydraulic conductivity of root xylem (in kg H2O·s-1·m-2).}
#'         \item{\code{VCleaf_kmax}: Maximum leaf hydraulic conductance.}
#'         \item{\code{VCleaf_c}, \code{VCleaf_d}: Parameters of the leaf vulnerability curve.}
#'         \item{\code{VCstem_kmax}: Maximum stem xylem conductance.}
#'         \item{\code{VCstem_c}, \code{VCstem_d}: Parameters of the stem xylem vulnerability curve.}
#'         \item{\code{VCroot_c}, \code{VCroot_d}: Parameters of the root xylem vulnerability curve.}
#'         \item{\code{Plant_kmax}: Maximum whole-plant conductance.}
#'       }
#'     }
#'     \item{\code{paramsWaterStorage}: A data frame with plant water storage parameters for each cohort:
#'       \itemize{
#'         \item{\code{LeafPI0}: Osmotic potential at full turgor of leaves (MPa).}
#'         \item{\code{LeafEPS}: Modulus of elasticity (capacity of the cell wall to resist changes in volume in response to changes in turgor) of leaves (MPa).}
#'         \item{\code{LeafAF}: Apoplastic fraction (proportion of water outside the living cells) in leaves.}
#'         \item{\code{Vleaf}: Storage water capacity in leaves, per leaf area (L/m2).}
#'         \item{\code{StemPI0}: Osmotic potential at full turgor of symplastic xylem tissue (MPa).}
#'         \item{\code{StemEPS}: Modulus of elasticity (capacity of the cell wall to resist changes in volume in response to changes in turgor) of symplastic xylem tissue (Mpa).}
#'         \item{\code{StemAF}: Apoplastic fraction (proportion of water outside the living cells) in stem xylem.}
#'         \item{\code{Vstem}: Storage water capacity in sapwood, per leaf area (L/m2).}
#'       }
#'     }
#'     \item{\code{internalPhenology} and \code{internalWater}: data frames to store internal state variables.}
#'     \item{\code{internalFCCS}: A data frame with fuel characteristics, according to \code{\link{fuel_FCCS}} (only if \code{fireHazardResults = TRUE}, in the control list).}
#'   }
#'   
#' Function \code{forest2growthInput} returns a list of class \code{growthInput} with the same elements as \code{spwbInput}, but with additional information. 
#' \itemize{
#' \item{Element \code{above} includes the following additional columns:
#'     \itemize{
#'       \item{\code{LA_live}: Live leaf area per individual (m2/ind).}
#'       \item{\code{LA_dead}: Dead leaf area per individual (m2/ind).}
#'       \item{\code{SA}: Live sapwood area per individual (cm2/ind).} 
#'   }
#'   }
#'   \item{\code{paramsGrowth}: A data frame with growth parameters for each cohort:
#'     \itemize{
#'       \item{\code{RERleaf}: Maintenance respiration rates (at 20ºC) for leaves (in g gluc·g dry-1·day-1).}
#'       \item{\code{RERsapwood}: Maintenance respiration rates (at 20ºC) for sapwood (in g gluc·g dry-1·day-1).}
#'       \item{\code{RERfineroot}: Maintenance respiration rates (at 20ºC) for fine roots (in g gluc·g dry-1·day-1).}
#'       \item{\code{CCleaf}: Leaf construction costs (in g gluc·g dry-1).}
#'       \item{\code{CCsapwood}: Sapwood construction costs (in g gluc·g dry-1).}
#'       \item{\code{CCfineroot}: Fine root construction costs (in g gluc·g dry-1).}
#'       \item{\code{RGRleafmax}: Maximum leaf relative growth rate (in m2·cm-2·day-1).}
#'       \item{\code{RGRsapwoodmax}: Maximum sapwood relative growth rate (in cm2·cm-2·day-1).}
#'       \item{\code{RGRfinerootmax}: Maximum fine root relative growth rate (in g dry·g dry-1·day-1).}
#'       \item{\code{SRsapwood}: Sapwood daily senescence rate (in day-1).}
#'       \item{\code{SRfineroot}: Fine root daily senescence rate (in day-1).}
#'       \item{\code{RSSG}: Minimum relative starch for sapwood growth (proportion).}
#'       \item{\code{fHDmin}: Minimum value of the height-to-diameter ratio (dimensionless).}
#'       \item{\code{fHDmax}: Maximum value of the height-to-diameter ratio (dimensionless).}
#'       \item{\code{WoodC}: Wood carbon content per dry weight (g C /g dry).}
#'     }
#'   }
#'   \item{\code{paramsMortalityRegeneration}: A data frame with mortality/regeneration parameters for each cohort:
#'     \itemize{
#'       \item{\code{MortalityBaselineRate}: Deterministic proportion or probability specifying the baseline reduction of cohort's density occurring in a year.}
#'       \item{\code{SurvivalModelStep}: Time step in years of the empirical survival model depending on stand basal area (e.g. 10).}
#'       \item{\code{SurvivalB0}: Intercept of the logistic baseline survival model depending on stand basal area.}
#'       \item{\code{SurvivalB1}: Slope of the logistic baseline survival model depending on stand basal area.}
#'       \item{\code{RecrTreeDensity}: Density of tree recruits from seeds.}
#'       \item{\code{IngrowthTreeDensity}: Density of trees reaching ingrowth DBH.}
#'       \item{\code{RecrTreeDBH}: DBH for tree recruits from seeds or resprouting (e.g. 1 cm).}
#'       \item{\code{IngrowthTreeDBH}: Ingrowth DBH for trees (e.g. 7.5 cm).}
#'     }
#'   }
#'   \item{\code{paramsAllometry}: A data frame with allometric parameters for each cohort:
#'     \itemize{
#'       \item{\code{Aash}: Regression coefficient relating the square of shrub height with shrub area.}
#'       \item{\code{Absh}, \code{Bbsh}: Allometric coefficients relating phytovolume with dry weight of shrub individuals.}
#'       \item{\code{Acr}, \code{B1cr}, \code{B2cr}, \code{B3cr}, \code{C1cr}, \code{C2cr}: Regression coefficients used to calculate crown ratio of trees.}
#'       \item{\code{Acw}, \code{Bcw}: Regression coefficients used to calculated crown width of trees.}
#'     }
#'   }
#'   \item {\code{internalAllocation}: A data frame with internal allocation variables for each cohort:
#'     \itemize{
#'       \item{\code{allocationTarget}: Value of the allocation target variable.}
#'       \item{\code{leafAreaTarget}: Target leaf area (m2) per individual.}
#'       \item{\code{sapwoodAreaTarget}: Target sapwood area (cm2) per individual.}
#'       \item{\code{fineRootBiomassTarget}: Target fine root biomass (g dry) per individual.}
#'       \item{\code{crownBudPercent}: Percentage of the crown with buds.}
#'     }
#'   }
#'   \item{\code{internalCarbon}: A data frame with the concentration (mol·gluc·l-1) of metabolic and storage carbon compartments for leaves and sapwood.}
#'   \item{\code{internalMortality}: A data frame to store the cumulative mortality (density for trees and cover for shrubs) predicted during the simulation,
#'   also distinguishing mortality due to starvation or dessication.}
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{resetInputs}}, \code{\link{spwb}}, \code{\link{soil}},  
#' \code{\link{forest}}, \code{\link{SpParamsMED}}, \code{\link{defaultSoilParams}}, \code{\link{plant_ID}}
#' 
#' @examples
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' # Aboveground parameters
#' forest2aboveground(exampleforestMED, SpParamsMED)
#' 
#' # Example of aboveground parameters taken from a forest
#' # described using LAI and crown ratio
#' data(exampleforestMED2)
#' forest2aboveground(exampleforestMED2, SpParamsMED)
#' 
#' # Initialize soil with default soil params
#' examplesoil <- soil(defaultSoilParams())
#' 
#' # Bewowground parameters (distribution of fine roots)
#' forest2belowground(exampleforestMED, examplesoil, SpParamsMED)
#' 
#' # Initialize control parameters using 'Granier' transpiration mode
#' control <- defaultControl("Granier")
#' 
#' # Prepare spwb input
#' forest2spwbInput(exampleforestMED, examplesoil, SpParamsMED, control)
#'                 
#' # Prepare input for 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' forest2spwbInput(exampleforestMED,examplesoil,SpParamsMED, control)
#' 
#' # Prepare input for 'Cochard' transpiration mode
#' control <- defaultControl("Cochard")
#' forest2spwbInput(exampleforestMED,examplesoil,SpParamsMED, control)
#' 
#' # Example of initialization from a forest 
#' # described using LAI and crown ratio
#' control <- defaultControl("Granier")
#' forest2spwbInput(exampleforestMED2, examplesoil, SpParamsMED, control)
#' 
#' @name modelInput
#' @aliases spwbInput growthInput
forest2spwbInput <- function(x, soil, SpParams, control) {
    .Call(`_medfate_forest2spwbInput`, x, soil, SpParams, control)
}

#' @rdname modelInput
forest2growthInput <- function(x, soil, SpParams, control) {
    .Call(`_medfate_forest2growthInput`, x, soil, SpParams, control)
}

#' Reset simulation inputs
#' 
#' Function \code{resetInputs()} allows resetting state variables in \code{x} to their defaults.
#' 
#' @param x An object of class \code{\link{spwbInput}} or \code{\link{growthInput}}.
#' 
#' @return Does not return any value. Instead, it modifies input object \code{x}.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{spwbInput}}, \code{\link{growthInput}}, \code{\link{spwb}}
#' 
resetInputs <- function(x) {
    invisible(.Call(`_medfate_resetInputs`, x))
}

.updateBelow <- function(x) {
    invisible(.Call(`_medfate_updateBelow`, x))
}

.multiplyInputParam <- function(x, paramType, paramName, cohort, f, message) {
    invisible(.Call(`_medfate_multiplyInputParam`, x, paramType, paramName, cohort, f, message))
}

.modifyInputParam <- function(x, paramType, paramName, cohort, newValue, message) {
    invisible(.Call(`_medfate_modifyInputParam`, x, paramType, paramName, cohort, newValue, message))
}

.checkSpeciesParameters <- function(SpParams, params) {
    invisible(.Call(`_medfate_checkSpeciesParameters`, SpParams, params))
}

.speciesNumericParameterFromSpIndex <- function(SP, SpParams, parName) {
    .Call(`_medfate_speciesNumericParameterFromIndex`, SP, SpParams, parName)
}

.speciesCharacterParameterFromSpIndex <- function(SP, SpParams, parName) {
    .Call(`_medfate_speciesCharacterParameterFromIndex`, SP, SpParams, parName)
}

#' @rdname species_values
species_characterParameter <- function(species, SpParams, parName) {
    .Call(`_medfate_speciesCharacterParameter`, species, SpParams, parName)
}

#' @rdname plant_values
plant_characterParameter <- function(x, SpParams, parName) {
    .Call(`_medfate_cohortCharacterParameter`, x, SpParams, parName)
}

#' @rdname species_values
species_parameter <- function(species, SpParams, parName, fillMissing = TRUE) {
    .Call(`_medfate_speciesNumericParameterWithImputation`, species, SpParams, parName, fillMissing)
}

#' @rdname plant_values
plant_parameter <- function(x, SpParams, parName, fillMissing = TRUE) {
    .Call(`_medfate_cohortNumericParameterWithImputation`, x, SpParams, parName, fillMissing)
}

.gdd <- function(DOY, Temp, Tbase = 5.0, cum = 0.0) {
    .Call(`_medfate_gdd`, DOY, Temp, Tbase, cum)
}

#' Leaf phenology
#'
#' Function \code{pheno_leafDevelopmentStatus} returns the expanded status (0 to 1) of leaves according to the growth degree days required to start bud burst and leaf unfolding, as dictated by a simple ecodormancy (one-phase) model (Chuine et al. 2013). 
#' Function \code{pheno_leafSenescenceStatus} returns the 0/1 senescence status of leaves according to the one-phase senescence model of Delpierre et al. (2009) on the basis of photoperiod and temperature.
#' Function \code{pheno_updateLeaves} updates the status of expanded leaves and dead leaves of object \code{x} given the photoperiod, temperature and wind of a given day. It applies the development model for 1 < doy < 180 and the senescence model for 181 > doy > 365.
#' 
#' @param Sgdd Degree days required for leaf budburst (in Celsius).
#' @param gdd Cumulative degree days (in Celsius)
#' @param unfoldingDD Degree-days for complete leaf unfolding after budburst has occurred.
#' 
#' @return Function \code{pheno_leafDevelopmentStatus} returns a vector of values between 0 and 1, 
#' whereas function \code{pheno_leafSenescenceStatus} returns a vector of 0 (senescent) and 1 (expanded) values. 
#' The other two functions do not return any value (see note).
#' 
#' @note Functions \code{pheno_updatePhenology} and \code{pheno_updateLeaves} modify the input object \code{x}. The first modifies phenological state and the second modifies the leaf area accordingly.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @references
#' Chuine, I., De Cortazar-Atauri, I.G., Kramer, K., \enc{Hänninen}{Hanninen}, H., 2013. Plant development models. Phenology: An Integrative Environmental Science. Springer, pp. 275–293.
#' 
#' Delpierre N, Dufrêne E, Soudani K et al (2009) Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in France. Agric For Meteorol 149:938–948. doi:10.1016/j.agrformet.2008.11.014
#' 
#' @seealso \code{\link{spwb}}, \code{\link{spwbInput}}
#' 
#' @name pheno_updateLeaves
pheno_leafDevelopmentStatus <- function(Sgdd, gdd, unfoldingDD = 300.0) {
    .Call(`_medfate_leafDevelopmentStatus`, Sgdd, gdd, unfoldingDD)
}

#' @param Ssen Threshold to start leaf senescence.
#' @param sen Cumulative senescence variable.
#' @rdname pheno_updateLeaves
pheno_leafSenescenceStatus <- function(Ssen, sen) {
    .Call(`_medfate_leafSenescenceStatus`, Ssen, sen)
}

#' @param x An object of class \code{\link{spwbInput}}.
#' @param doy Day of the year.
#' @param photoperiod Day length (in hours).
#' @param tmean Average day temperature (in Celsius).
#' @rdname pheno_updateLeaves
pheno_updatePhenology <- function(x, doy, photoperiod, tmean) {
    invisible(.Call(`_medfate_updatePhenology`, x, doy, photoperiod, tmean))
}

#' @param wind Average day wind speed (in m/s).
#' @param fromGrowthModel Boolean flag to indicate that routine is called from \code{\link{growth}} simulation function.
#' @rdname pheno_updateLeaves
pheno_updateLeaves <- function(x, wind, fromGrowthModel) {
    invisible(.Call(`_medfate_updateLeaves`, x, wind, fromGrowthModel))
}

#' Photosynthesis submodel functions
#' 
#' Set of functions used in the calculation of photosynthesis
#' 
#' @param Tleaf Leaf temperature (in ºC).
#' @param Oi Oxigen concentration (mmol*mol-1).
#' @param Vmax298,Vmax298SL,Vmax298SH Maximum Rubisco carboxylation rate per leaf area at 298ºK (i.e. 25 ºC) (micromol*s-1*m-2) (for each canopy layer in the case of \code{photo_multilayerPhotosynthesisFunction}). 'SH' stands for shade leaves, whereas 'SL' stands for sunlit leaves.
#' @param Jmax298,Jmax298SL,Jmax298SH Maximum electron transport rate per leaf area at 298ºK (i.e. 25 ºC) (micromol*s-1*m-2) (for each canopy layer in the case of \code{photo_multilayerPhotosynthesisFunction}). 'SH' stands for shade leaves, whereas 'SL' stands for sunlit leaves.
#' @param Q Active photon flux density (micromol * s-1 * m-2).
#' @param Ci CO2 internal concentration (micromol * mol-1).
#' @param GT CO2 saturation point corrected by temperature (micromol * mol-1).
#' @param Jmax Maximum electron transport rate per leaf area (micromol*s-1*m-2).
#' @param Km Km = Kc*(1.0+(Oi/Ko)) - Michaelis-Menten term corrected by temperature (in micromol * mol-1).
#' @param Vmax Maximum Rubisco carboxylation rate per leaf area (micromol*s-1*m-2).
#' @param Catm CO2 air concentration (micromol * mol-1).
#' @param Gc CO2 leaf (stomatal) conductance (mol * s-1 * m-2).
#' @param E Transpiration flow rate per leaf area (mmol*s-1*m-2).
#' @param psiLeaf Leaf water potential (MPa).
#' @param Patm Atmospheric air pressure (in kPa).
#' @param Tair Air temperature (in ºC).
#' @param vpa Vapour pressure deficit (in kPa).
#' @param u Wind speed above the leaf boundary (in m/s) (for each canopy layer in the case of \code{photo_multilayerPhotosynthesisFunction}).
#' @param absRad Absorbed long- and short-wave radiation (in W*m^-2).
#' @param SWRabs Absorbed short-wave radiation (in W·m-2).
#' @param LWRnet Net long-wave radiation balance (in W·m-2).
#' @param leafWidth Leaf width (in cm).
#' @param refLeafArea Leaf reference area.
#' @param verbose Boolean flag to indicate console output.
#' @param SLarea,SHarea Leaf area index of sunlit/shade leaves (for each canopy layer in the case of \code{photo_multilayerPhotosynthesisFunction}).
#' @param absRadSL,absRadSH Instantaneous absorbed radiation (W·m-2) per unit of sunlit/shade leaf area (for each canopy layer in the case of \code{photo_multilayerPhotosynthesisFunction}).
#' @param QSL,QSH Active photon flux density (micromol * s-1 * m-2) per unit of sunlit/shade leaf area (for each canopy layer in the case of \code{photo_multilayerPhotosynthesisFunction}).
#' 
#' @details Details of the photosynthesis submodel are given in the medfate book
#' 
#' @return
#' Values returned for each function are:
#' \itemize{
#'   \item{\code{photo_GammaTemp}: CO2 compensation concentration (micromol * mol-1).}
#'   \item{\code{photo_KmTemp}: Michaelis-Menten coefficients of Rubisco for Carbon (micromol * mol-1) and Oxigen (mmol * mol-1).}
#'   \item{\code{photo_VmaxTemp}: Temperature correction of Vmax298.}
#'   \item{\code{photo_JmaxTemp}: Temperature correction of Jmax298.}
#'   \item{\code{photo_electronLimitedPhotosynthesis}: Electron-limited photosynthesis (micromol*s-1*m-2) following Farquhar et al. (1980).}
#'   \item{\code{photo_rubiscoLimitedPhotosynthesis}: Rubisco-limited photosynthesis (micromol*s-1*m-2) following Farquhar et al. (1980).}
#'   \item{\code{photo_photosynthesis}: Calculates gross photosynthesis (micromol*s-1*m-2) following (Farquhar et al. (1980) and Collatz et al (1991).}
#'   \item{\code{photo_leafPhotosynthesisFunction}: Returns a data frame with the following columns:
#'     \itemize{
#'       \item{\code{LeafTemperature}: Leaf temperature (ºC).}
#'       \item{\code{LeafVPD}: Leaf vapor pressure deficit (kPa).}
#'       \item{\code{LeafCi}: Internal CO2 concentration (micromol * mol-1).}
#'       \item{\code{Gsw}: Leaf stomatal conductance to water vapor (mol * s-1 * m-2).}
#'       \item{\code{GrossPhotosynthesis}: Gross photosynthesis (micromol*s-1*m-2).}
#'       \item{\code{NetPhotosynthesis}: Net photosynthesis, after discounting autotrophic respiration (micromol*s-1*m-2).}
#'     }
#'   }
#'   \item{\code{photo_sunshadePhotosynthesisFunction}: Returns a data frame with the following columns:
#'     \itemize{
#'       \item{\code{GrossPhotosynthesis}: Gross photosynthesis (micromol*s-1*m-2).}
#'       \item{\code{NetPhotosynthesis}: Net photosynthesis, after discounting autotrophic respiration (micromol*s-1*m-2).}
#'       \item{\code{LeafCiSL}: Sunlit leaf internal CO2 concentration (micromol * mol-1).}
#'       \item{\code{LeafCiSH}: Shade leaf internal CO2 concentration (micromol * mol-1).}
#'       \item{\code{LeafTempSL}: Sunlit leaf temperature (ºC).}
#'       \item{\code{LeafTempSH}: Shade leaf temperature (ºC).}
#'       \item{\code{LeafVPDSL}: Sunlit leaf vapor pressure deficit (kPa).}
#'       \item{\code{LeafVPDSH}: Shade leaf vapor pressure deficit (kPa).}
#'     }
#'   }
#'   \item{\code{photo_multilayerPhotosynthesisFunction}: Return a data frame with the following columns:
#'     \itemize{
#'       \item{\code{GrossPhotosynthesis}: Gross photosynthesis (micromol*s-1*m-2).}
#'       \item{\code{NetPhotosynthesis}: Net photosynthesis, after discounting autotrophic respiration (micromol*s-1*m-2).}
#'     }
#'   }
#' }
#' 
#' @references
#' Bernacchi, C. J., E. L. Singsaas, C. Pimentel, A. R. Portis, and S. P. Long. 2001. Improved temperature response functions for models of Rubisco-limited photosynthesis. Plant, Cell and Environment 24:253–259.
#' 
#' Collatz, G. J., J. T. Ball, C. Grivet, and J. A. Berry. 1991. Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer. Agricultural and Forest Meteorology 54:107–136.
#' 
#' Farquhar, G. D., S. von Caemmerer, and J. A. Berry. 1980. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149:78–90.
#' 
#' Leuning, R. 2002. Temperature dependence of two parameters in a photosynthesis model. Plant, Cell and Environment 25:1205–1210.
#' 
#' Sperry, J. S., M. D. Venturas, W. R. L. Anderegg, M. Mencuccini, D. S. Mackay, Y. Wang, and D. M. Love. 2016. Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost. Plant Cell and Environment.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{hydraulics_supplyFunctionNetwork}}, \code{\link{biophysics_leafTemperature}}, \code{\link{spwb}}
#' 
#' @name photo
photo_GammaTemp <- function(Tleaf) {
    .Call(`_medfate_gammaTemp`, Tleaf)
}

#' @rdname photo
photo_KmTemp <- function(Tleaf, Oi = 209.0) {
    .Call(`_medfate_KmTemp`, Tleaf, Oi)
}

#' @rdname photo
photo_VmaxTemp <- function(Vmax298, Tleaf) {
    .Call(`_medfate_VmaxTemp`, Vmax298, Tleaf)
}

#' @rdname photo
photo_JmaxTemp <- function(Jmax298, Tleaf) {
    .Call(`_medfate_JmaxTemp`, Jmax298, Tleaf)
}

#' @rdname photo
photo_electronLimitedPhotosynthesis <- function(Q, Ci, GT, Jmax) {
    .Call(`_medfate_electronLimitedPhotosynthesis`, Q, Ci, GT, Jmax)
}

#' @rdname photo
photo_rubiscoLimitedPhotosynthesis <- function(Ci, GT, Km, Vmax) {
    .Call(`_medfate_rubiscoLimitedPhotosynthesis`, Ci, GT, Km, Vmax)
}

#' @rdname photo
photo_photosynthesis <- function(Q, Catm, Gc, Tleaf, Vmax298, Jmax298, verbose = FALSE) {
    .Call(`_medfate_leafphotosynthesis`, Q, Catm, Gc, Tleaf, Vmax298, Jmax298, verbose)
}

#' @rdname photo
photo_leafPhotosynthesisFunction <- function(E, psiLeaf, Catm, Patm, Tair, vpa, u, absRad, Q, Vmax298, Jmax298, leafWidth = 1.0, refLeafArea = 1.0, verbose = FALSE) {
    .Call(`_medfate_leafPhotosynthesisFunction`, E, psiLeaf, Catm, Patm, Tair, vpa, u, absRad, Q, Vmax298, Jmax298, leafWidth, refLeafArea, verbose)
}

#' @rdname photo
photo_leafPhotosynthesisFunction2 <- function(E, psiLeaf, Catm, Patm, Tair, vpa, u, SWRabs, LWRnet, Q, Vmax298, Jmax298, leafWidth = 1.0, refLeafArea = 1.0, verbose = FALSE) {
    .Call(`_medfate_leafPhotosynthesisFunction2`, E, psiLeaf, Catm, Patm, Tair, vpa, u, SWRabs, LWRnet, Q, Vmax298, Jmax298, leafWidth, refLeafArea, verbose)
}

#' @rdname photo
photo_sunshadePhotosynthesisFunction <- function(E, psiLeaf, Catm, Patm, Tair, vpa, SLarea, SHarea, u, absRadSL, absRadSH, QSL, QSH, Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH, leafWidth = 1.0, verbose = FALSE) {
    .Call(`_medfate_sunshadePhotosynthesisFunction`, E, psiLeaf, Catm, Patm, Tair, vpa, SLarea, SHarea, u, absRadSL, absRadSH, QSL, QSH, Vmax298SL, Vmax298SH, Jmax298SL, Jmax298SH, leafWidth, verbose)
}

#' @rdname photo
photo_multilayerPhotosynthesisFunction <- function(E, psiLeaf, Catm, Patm, Tair, vpa, SLarea, SHarea, u, absRadSL, absRadSH, QSL, QSH, Vmax298, Jmax298, leafWidth = 1.0, verbose = FALSE) {
    .Call(`_medfate_multilayerPhotosynthesisFunction`, E, psiLeaf, Catm, Patm, Tair, vpa, SLarea, SHarea, u, absRadSL, absRadSH, QSL, QSH, Vmax298, Jmax298, leafWidth, verbose)
}

#' Root functions
#' 
#' Functions to calculate properties of fine/coarse roots within the soil, given root system parameters and soil layer definition.
#' 
#' @param Z50 A vector of depths (in mm) corresponding to 50\% of roots.
#' @param Z95 A vector of depths (in mm) corresponding to 95\% of roots.
#' @param Zcone A vector of depths (in mm) corresponding to the root cone tip.
#' @param d The width (in mm) corresponding to each soil layer.
#' @param v Vector of proportions of fine roots in each soil layer.
#' @param depthWidthRatio Ratio between radius of the soil layer with the largest radius and maximum rooting depth.
#' @param rfc Percentage of rock fragment content (volume basis) for each layer.
#' @param Kmax_rootxylem Sapwood-specific hydraulic conductivity of root xylem (in kg H2O·s-1·m-1·MPa-1).
#' @param VCroot_kmax Root xylem maximum conductance per leaf area (mmol·m-2·s-1·MPa-1). 
#' @param Al2As Leaf area to sapwood area ratio (in m2·m-2).
#' @param specificRootLength Specific fine root length (length of fine roots over weight).
#' @param rootTissueDensity Fine root tissue density (weight over volume at turgidity).
#' @param Ksoil Soil saturated conductivity (mmol·m-1·s-1·MPa-1).
#' @param krhizo Rhizosphere maximum conductance per leaf area (mmol·m-2·s-1·MPa-1).
#' @param lai Leaf area index.
#' @param rootLengthDensity Fine root length density (length of fine roots over soil volume; cm/cm3)
#' @param fineRootBiomass Biomass of fine roots (g).
#' @param V Matrix of proportions of fine roots (cohorts x soil layers).
#' @param VolInd Volume of soil (in m3) occupied by coarse roots per individual. 
#' @param N Density of individuals per hectare.
#' @param poolProportions Division of the stand area among plant cohorts (proportions).
#' 
#' @details
#' \itemize{
#'   \item{\code{root_conicDistribution()} assumes a (vertical) conic distribution of fine roots, whereas \code{root_ldrDistribution()} distributes fine roots according to the linear dose response model of Schenck & Jackson (2002). Return a matrix of fine root proportions in each layer with as many rows as elements in \code{Z} (or \code{Z50}) and as many columns as soil layers.}
#'   \item{\code{root_coarseRootLengths()} and \code{root_coarseRootLengthsFromVolume()} estimate the length of coarse roots (mm) for each soil layer, including axial and radial lengths.}
#'   \item{\code{root_coarseRootSoilVolume} estimates the soil volume (m3) occupied by coarse roots of an individual.}
#'   \item{\code{root_coarseRootSoilVolumeFromConductance} estimates the soil volume (m3) occupied by coarse roots of an individual from root xylem conductance.}
#'   \item{\code{root_fineRootHalfDistance()} calculates the half distance (cm) between neighbouring fine roots.}
#'   \item{\code{root_fineRootRadius()} calculates the radius of fine roots (cm).}
#'   \item{\code{root_fineRootAreaIndex()} estimates the fine root area index for a given soil conductivity and maximum rhizosphere conductance.}
#'   \item{\code{root_fineRootBiomass()} estimates the biomass of fine roots (g dry/individual) for a given soil conductivity and maximum rhizosphere conductance.}
#'   \item{\code{root_rhizosphereMaximumConductance()} is the inverse of the preceeding function, i.e. it estimates rhizosphere conductance from soil conductivity and fine root biomass.}
#'   \item{\code{root_fineRootSoilVolume()} calculates the soil volume (m3) occupied with fine roots.}
#'   \item{\code{root_specificRootSurfaceArea()} returns the specific fine root area (cm2/g).}
#'   \item{\code{root_individualRootedGroundArea()} calculates the area (m2) covered by roots of an individual, for each soil layer.}
#'   \item{\code{root_horizontalProportions()} calculates the (horizontal) proportion of roots of each cohort in the water pool corresponding to itself and that of other cohorts, for each soil layer. Returns a list (with as many elements as cohorts) with each element being a matrix.}
#'   }
#' 
#' @return See details.
#' 
#' @references
#' Schenk, H., Jackson, R., 2002. The global biogeography of roots. Ecol. Monogr. 72, 311–328.
#' 
#' Sperry, J. S., Y. Wang, B. T. Wolfe, D. S. Mackay, W. R. L. Anderegg, N. G. Mcdowell, and W. T. Pockman. 2016. Pragmatic hydraulic theory predicts stomatal responses to climatic water deficits. New Phytologist 212, 577–589.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#'  \code{\link{spwb}},  \code{\link{forest2spwbInput}}, \code{\link{soil}}
#'
#' @examples
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' ntree = nrow(exampleforestMED$treeData)
#' 
#' #Initialize soil with default soil params
#' S = soil(defaultSoilParams())
#' 
#' #Calculate conic root system for trees
#' V1 = root_conicDistribution(Z=rep(2000,ntree), S$dVec)            
#' print(V1)
#'      
#' #Calculate LDR root system for trees (Schenck & Jackson 2002)
#' V2 = root_ldrDistribution(Z50 = rep(200,ntree), 
#'                           Z95 = rep(1000,ntree), S$dVec)
#' print(V2)     
#' 
#' @name root
root_conicDistribution <- function(Zcone, d) {
    .Call(`_medfate_conicDistribution`, Zcone, d)
}

#' @rdname root
root_ldrDistribution <- function(Z50, Z95, d) {
    .Call(`_medfate_ldrDistribution`, Z50, Z95, d)
}

.rootDistribution <- function(z, x) {
    .Call(`_medfate_rootDistribution`, z, x)
}

#' @rdname root
root_individualRootedGroundArea <- function(VolInd, V, d, rfc) {
    .Call(`_medfate_individualRootedGroundArea`, VolInd, V, d, rfc)
}

#' @rdname root
root_specificRootSurfaceArea <- function(specificRootLength, rootTissueDensity) {
    .Call(`_medfate_specificRootSurfaceArea`, specificRootLength, rootTissueDensity)
}

#' @rdname root
root_fineRootRadius <- function(specificRootLength, rootTissueDensity) {
    .Call(`_medfate_fineRootRadius`, specificRootLength, rootTissueDensity)
}

#' @rdname root
root_fineRootHalfDistance <- function(rootLengthDensity) {
    .Call(`_medfate_fineRootHalfDistance`, rootLengthDensity)
}

#' @rdname root
root_fineRootAreaIndex <- function(Ksoil, krhizo, lai, specificRootLength, rootTissueDensity, rootLengthDensity) {
    .Call(`_medfate_fineRootAreaIndex`, Ksoil, krhizo, lai, specificRootLength, rootTissueDensity, rootLengthDensity)
}

#' @rdname root
root_fineRootBiomass <- function(Ksoil, krhizo, lai, N, specificRootLength, rootTissueDensity, rootLengthDensity) {
    .Call(`_medfate_fineRootBiomassPerIndividual`, Ksoil, krhizo, lai, N, specificRootLength, rootTissueDensity, rootLengthDensity)
}

#' @rdname root
root_rhizosphereMaximumConductance <- function(Ksoil, fineRootBiomass, lai, N, specificRootLength, rootTissueDensity, rootLengthDensity) {
    .Call(`_medfate_rhizosphereMaximumConductance`, Ksoil, fineRootBiomass, lai, N, specificRootLength, rootTissueDensity, rootLengthDensity)
}

#' @rdname root
root_fineRootSoilVolume <- function(fineRootBiomass, specificRootLength, rootLengthDensity) {
    .Call(`_medfate_fineRootSoilVolume`, fineRootBiomass, specificRootLength, rootLengthDensity)
}

#' @rdname root
root_coarseRootSoilVolumeFromConductance <- function(Kmax_rootxylem, VCroot_kmax, Al2As, v, d, rfc) {
    .Call(`_medfate_coarseRootSoilVolumeFromConductance`, Kmax_rootxylem, VCroot_kmax, Al2As, v, d, rfc)
}

#' @rdname root
root_coarseRootLengthsFromVolume <- function(VolInd, v, d, rfc) {
    .Call(`_medfate_coarseRootLengthsFromVolume`, VolInd, v, d, rfc)
}

#' @rdname root
root_coarseRootLengths <- function(v, d, depthWidthRatio = 1.0) {
    .Call(`_medfate_coarseRootLengths`, v, d, depthWidthRatio)
}

#' @rdname root
root_coarseRootSoilVolume <- function(v, d, depthWidthRatio = 1.0) {
    .Call(`_medfate_coarseRootSoilVolume`, v, d, depthWidthRatio)
}

#' @rdname root
root_horizontalProportions <- function(poolProportions, VolInd, N, V, d, rfc) {
    .Call(`_medfate_horizontalProportions`, poolProportions, VolInd, N, V, d, rfc)
}

#' Soil texture and hydraulics
#' 
#' Low-level functions relating soil texture with soil hydraulics and soil water content.
#'
#' @param clay Percentage of clay (in percent weight).
#' @param sand Percentage of sand (in percent weight).
#' @param n,alpha,theta_res,theta_sat Parameters of the Van Genuchten-Mualem model (m = 1 - 1/n).
#' @param psi Water potential (in MPa).
#' @param theta Relative water content (in percent volume).
#' @param om Percentage of organic matter (optional, in percent weight).
#' @param mmol Boolean flag to indicate that saturated conductivity units should be returned in mmol/m/s/MPa. If \code{mmol = FALSE} then units are cm/day.
#' @param bd Bulk density (in g/cm3).
#' @param topsoil A boolean flag to indicate topsoil layer.
#' @param soilType A string indicating the soil type.
#' @param soil Soil object (returned by function \code{\link{soil}}).
#' @param model Either 'SX' or 'VG' for Saxton's or Van Genuchten's water retention models; or 'both' to plot both retention models.
#' @param minPsi Minimum water potential (in MPa) to calculate the amount of extractable water.
#' @param pWeight Percentage of corresponding to rocks, in weight.
#' @param bulkDensity Bulk density of the soil fraction (g/cm3).
#' @param rockDensity Rock density (g/cm3).
#' 
#' @details
#' \itemize{
#' \item{\code{soil_psi2thetaSX()} and \code{soil_theta2psiSX()} calculate water potentials (MPa) and water contents (theta) using texture data the formulae of Saxton et al. (1986) or Saxton & Rawls (2006) depending on whether organic matter is available.}
#' \item{\code{soil_psi2thetaVG()} and \code{soil_theta2psiVG()} to the same calculations as before, but using the Van Genuchten - Mualem equations (\enc{Wösten}{Wosten} & van Genuchten 1988). }
#' \item{\code{soil_saturatedConductivitySX()} returns the saturated conductivity of the soil (in cm/day or mmol/m/s/MPa), estimated from formulae of Saxton et al. (1986) or Saxton & Rawls (2006) depending on whether organic matter is available.}
#' \item{\code{soil_unsaturatedConductivitySX()} returns the unsaturated conductivity of the soil (in cm/day or mmol/m/s/MPa), estimated from formulae of Saxton et al. (1986) or Saxton & Rawls (2006) depending on whether organic matter is available.}
#' \item{\code{soil_USDAType()} returns the USDA type (a string) for a given texture.}
#' \item{\code{soil_vanGenuchtenParamsCarsel()} gives parameters for van Genuchten-Mualem equations (alpha, n, theta_res and theta_sat, where alpha is in MPa-1) for a given texture type (Leij et al. 1996) }
#' \item{\code{soil_vanGenuchtenParamsToth()} gives parameters for van Genuchten-Mualem equations (alpha, n, theta_res and theta_sat, where alpha is in MPa-1) for a given texture, organic matter and bulk density (Toth et al. 2015).}
#' \item{\code{soil_psi()} returns the water potential (MPa) of each soil layer, according to its water retention model.}
#' \item{\code{soil_theta()} returns the moisture content (as percent of soil volume) of each soil layer, according to its water retention model.}
#' \item{\code{soil_water()} returns the water volume (mm) of each soil layer, according to its water retention model.}
#' \item{\code{soil_conductivity()} returns the conductivity of each soil layer (mmol/m/s/MPa), according the Saxton model.}
#' \item{\code{soil_waterExtractable()} returns the water volume (mm) extractable from the soil according to its water retention curves and up to a given soil water potential.}
#' \item{\code{soil_waterFC()} and \code{soil_thetaFC()} calculate the water volume (in mm) and moisture content (as percent of soil volume) of each soil layer at field capacity, respectively.}
#' \item{\code{soil_waterWP()} and \code{soil_thetaWP()} calculate the water volume (in mm) and moisture content (as percent of soil volume) of each soil layer at wilting point (-1.5 MPa), respectively. }
#' \item{\code{soil_waterSAT()}, \code{soil_thetaSATSX()} and \code{soil_thetaSAT()} calculate the saturated water volume (in mm) and moisture content (as percent of soil volume) of each soil layer.}
#' \item{\code{soil_waterTableDepth()} returns water table depth in mm from surface.}
#' \item{\code{soil_rockWeight2Volume()} transforms rock percentage from weight to volume basis.}
#' \item{\code{soil_retentionCurvePlot()} allows ploting the water retention curve of a given soil layer.}
#' }
#' 
#' @return Depends on the function (see details).
#' 
#' @references
#' Leij, F.J., Alves, W.J., Genuchten, M.T. Van, Williams, J.R., 1996. The UNSODA Unsaturated Soil Hydraulic Database User’s Manual Version 1.0.
#' 
#' Saxton, K.E., Rawls, W.J., Romberger, J.S., Papendick, R.I., 1986. Estimating generalized soil-water characteristics from texture. Soil Sci. Soc. Am. J. 50, 1031–1036.
#' 
#' Saxton, K.E., Rawls, W.J., 2006. Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil Sci. Soc. Am. J. 70, 1569. doi:10.2136/sssaj2005.0117
#' 
#' \enc{Wösten}{Wosten}, J.H.M., & van Genuchten, M.T. 1988. Using texture and other soil properties to predict the unsaturated soil hydraulic functions. Soil Science Society of America Journal 52: 1762–1770.
#' 
#' \enc{Tóth}{Toth}, B., Weynants, M., Nemes, A., \enc{Makó}{Mako}, A., Bilas, G., and \enc{Tóth}{Toth}, G. 2015. New generation of hydraulic pedotransfer functions for Europe. European Journal of Soil Science 66: 226–238.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso  \code{\link{soil}}
#' 
#' @examples
#' #Determine USDA soil texture type
#' type = soil_USDAType(clay=40, sand=10)
#' type
#' 
#' #Van Genuchten's params (bulk density = 1.3 g/cm)
#' vg = soil_vanGenuchtenParamsToth(40,10,1,1.3,TRUE)
#' vg
#' 
#' # Initialize soil object with default params
#' s = soil(defaultSoilParams())
#' 
#' # Plot Saxton's and Van Genuchten's water retention curves
#' soil_retentionCurvePlot(s, model="both")
#' 
#' @name soil_texture
soil_saturatedConductivitySX <- function(clay, sand, om = NA_real_, mmol = TRUE) {
    .Call(`_medfate_saturatedConductivitySaxton`, clay, sand, om, mmol)
}

#' @rdname soil_texture
soil_unsaturatedConductivitySX <- function(theta, clay, sand, om = NA_real_, mmol = TRUE) {
    .Call(`_medfate_unsaturatedConductivitySaxton`, theta, clay, sand, om, mmol)
}

#' @rdname soil_texture
soil_thetaSATSX <- function(clay, sand, om = NA_real_) {
    .Call(`_medfate_thetaSATSaxton`, clay, sand, om)
}

#' @rdname soil_texture
soil_theta2psiSX <- function(clay, sand, theta, om = NA_real_) {
    .Call(`_medfate_theta2psiSaxton`, clay, sand, theta, om)
}

#' @rdname soil_texture
soil_psi2thetaSX <- function(clay, sand, psi, om = NA_real_) {
    .Call(`_medfate_psi2thetaSaxton`, clay, sand, psi, om)
}

#' @rdname soil_texture
soil_psi2thetaVG <- function(n, alpha, theta_res, theta_sat, psi) {
    .Call(`_medfate_psi2thetaVanGenuchten`, n, alpha, theta_res, theta_sat, psi)
}

#' @rdname soil_texture
soil_theta2psiVG <- function(n, alpha, theta_res, theta_sat, theta) {
    .Call(`_medfate_theta2psiVanGenuchten`, n, alpha, theta_res, theta_sat, theta)
}

#' @rdname soil_texture
soil_USDAType <- function(clay, sand) {
    .Call(`_medfate_USDAType`, clay, sand)
}

#' @rdname soil_texture
soil_thetaFC <- function(soil, model = "SX") {
    .Call(`_medfate_thetaFC`, soil, model)
}

#' @rdname soil_texture
soil_thetaWP <- function(soil, model = "SX") {
    .Call(`_medfate_thetaWP`, soil, model)
}

#' @rdname soil_texture
soil_thetaSAT <- function(soil, model = "SX") {
    .Call(`_medfate_thetaSAT`, soil, model)
}

#' @rdname soil_texture
soil_waterFC <- function(soil, model = "SX") {
    .Call(`_medfate_waterFC`, soil, model)
}

#' @rdname soil_texture
soil_waterSAT <- function(soil, model = "SX") {
    .Call(`_medfate_waterSAT`, soil, model)
}

#' @rdname soil_texture
soil_waterWP <- function(soil, model = "SX") {
    .Call(`_medfate_waterWP`, soil, model)
}

#' @rdname soil_texture
soil_waterExtractable <- function(soil, model = "SX", minPsi = -5.0) {
    .Call(`_medfate_waterExtractable`, soil, model, minPsi)
}

#' @rdname soil_texture
soil_theta <- function(soil, model = "SX") {
    .Call(`_medfate_theta`, soil, model)
}

#' @rdname soil_texture
soil_water <- function(soil, model = "SX") {
    .Call(`_medfate_water`, soil, model)
}

#' @rdname soil_texture
soil_rockWeight2Volume <- function(pWeight, bulkDensity, rockDensity = 2.3) {
    .Call(`_medfate_rockWeight2Volume`, pWeight, bulkDensity, rockDensity)
}

#' @rdname soil_texture
soil_psi <- function(soil, model = "SX") {
    .Call(`_medfate_psi`, soil, model)
}

#' @rdname soil_texture
soil_conductivity <- function(soil) {
    .Call(`_medfate_conductivity`, soil)
}

#' @rdname soil_texture
soil_waterTableDepth <- function(soil, model = "SX") {
    .Call(`_medfate_waterTableDepth`, soil, model)
}

#' @rdname soil_texture
soil_vanGenuchtenParamsCarsel <- function(soilType) {
    .Call(`_medfate_vanGenuchtenParamsCarsel`, soilType)
}

#' @rdname soil_texture
soil_vanGenuchtenParamsToth <- function(clay, sand, om, bd, topsoil) {
    .Call(`_medfate_vanGenuchtenParamsToth`, clay, sand, om, bd, topsoil)
}

#' Soil initialization
#'
#' Initializes soil parameters and state variables for its use in simulations.
#' 
#' @param SoilParams A data frame of soil parameters (see an example in \code{\link{defaultSoilParams}}).
#' @param VG_PTF Pedotransfer functions to obtain parameters for the van Genuchten-Mualem equations. Either \code{"Carsel"} (Carsel and Parrish 1988) or \code{"Toth"} (Toth et al. 2015).
#' @param W A numerical vector with the initial relative water content of each soil layer.
#' @param SWE Initial snow water equivalent of the snow pack on the soil surface (mm).
#' 
#' @details 
#' Function \code{print} prompts a description of soil characteristics and state variables (water content and temperature) 
#' according to a water retention curve (either Saxton's or Van Genuchten's). 
#' Volume at field capacity is calculated assuming a soil water potential equal to -0.033 MPa. 
#' Parameter \code{Temp} is initialized as missing for all soil layers. 
#' 
#' @return
#' Function \code{soil} returns a list of class \code{soil} with the following elements:
#' \itemize{
#'   \item{\code{SoilDepth}: Soil depth (in mm).}
#'   \item{\code{W}: State variable with relative water content of each layer (in as proportion relative to FC).}
#'   \item{\code{Temp}: State variable with temperature (in ºC) of each layer.}
#'   \item{\code{Ksoil}: Kappa parameter for bare soil evaporation.}
#'   \item{\code{Gsoil}: Gamma parameter for bare soil evaporation (see \code{\link{hydrology_soilEvaporationAmount}}).}
#'   \item{\code{dVec}: Width of soil layers (in mm).}
#'   \item{\code{sand}: Sand percentage for each layer (in percent volume).}
#'   \item{\code{clay}: Clay percentage for each layer (in percent volume).}
#'   \item{\code{om}: Organic matter percentage for each layer (in percent volume).}
#'   \item{\code{VG_alpha}, \code{VG_n}, \code{VG_theta_res}, \code{VG_theta_sat}: Parameters for van Genuchten's pedotransfer functions, for each layer, corresponding to the USDA texture type.}
#'   \item{\code{Ksat}: Saturated soil conductivity for each layer (estimated using function \code{\link{soil_saturatedConductivitySX}}.}
#'   \item{\code{macro}: Macroporosity for each layer (estimated using Stolf et al. 2011).}
#'   \item{\code{rfc}: Percentage of rock fragment content for each layer.}
#'   \item{\code{Kdrain}: Saturated vertical hydraulic conductivity (mm/day) (i.e. how easy is deep drainage towards groundwater). Function \code{soil} estimates it as a function of soil saturated hydraulic conductivity, but should be parametrized as a function of bedrock material. }
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @references
#' Carsel, R.F., and Parrish, R.S. 1988. Developing joint probability distributions of soil water retention characteristics. Water Resources Research 24: 755–769.
#' 
#' \enc{Tóth}{Toth}, B., Weynants, M., Nemes, A., \enc{Makó}{Mako}, A., Bilas, G., and \enc{Tóth}{Toth}, G. 2015. New generation of hydraulic pedotransfer functions for Europe. European Journal of Soil Science 66: 226–238.
#' 
#' Stolf, R., Thurler, A., Oliveira, O., Bacchi, S., Reichardt, K., 2011. Method to estimate soil macroporosity and microporosity based on sand content and bulk density. Rev. Bras. Ciencias do Solo 35, 447–459.
#' 
#' @seealso   \code{\link{soil_psi2thetaSX}}, \code{\link{soil_psi2thetaVG}}, \code{\link{spwb}}, \code{\link{defaultSoilParams}}
#' 
#' @examples
#' # Initializes soil
#' s = soil(defaultSoilParams())
#' 
#' # Prints soil characteristics according to Saxton's water retention curve
#' print(s, model="SX")
#' 
#' # Prints soil characteristics according to Van Genuchten's water retention curve
#' print(s, model="VG")
#' 
#' @name soil
soil <- function(SoilParams, VG_PTF = "Toth", W = as.numeric( c(1.0)), SWE = 0.0) {
    .Call(`_medfate_soil`, SoilParams, VG_PTF, W, SWE)
}

.modifySoilLayerParam <- function(soil, paramName, layer, newValue, VG_PTF = "Toth") {
    invisible(.Call(`_medfate_modifySoilLayerParam`, soil, paramName, layer, newValue, VG_PTF))
}

#' Soil thermodynamic functions
#' 
#' Functions \code{soil_thermalConductivity} and \code{soil_thermalCapacity} calculate thermal conductivity and thermal capacity 
#' for each soil layer, given its texture and water content. Functions \code{soil_temperatureGradient} and \code{soil_temperatureChange} 
#' are used to calculate soil temperature gradients (in ºC/m) and temporal temperature change (in ºC/s) 
#' given soil layer texture and water content (and possibly including heat flux from above).
#' 
#' @param soil Soil object (returned by function \code{\link{soil}}).
#' @param model Either 'SX' or 'VG' for Saxton's or Van Genuchten's pedotransfer models.
#' @param dVec Width of soil layers (in mm).
#' @param Temp Temperature (in ºC) for each soil layer.
#' @param clay Percentage of clay (in percent weight) for each layer.
#' @param sand Percentage of sand (in percent weight) for each layer.
#' @param W Soil moisture (in percent of field capacity) for each layer.
#' @param Theta_FC Relative water content (in percent volume) at field capacity for each layer.
#' @param Gdown Downward heat flux from canopy to soil (in W·m-2).
#' 
#' @return 
#' Function \code{soil_thermalConductivity} returns a vector with values of thermal conductivity (W/m/ºK) for each soil layer. 
#' 
#' Function \code{soil_thermalCapacity} returns a vector with values of heat storage capacity (J/m3/ºK) for each soil layer. 
#' 
#' Function \code{soil_temperatureGradient} returns a vector with values of temperature gradient between consecutive soil layers. 
#' 
#' Function \code{soil_temperatureChange} returns a vector with values of instantaneous temperature change (ºC/s) for each soil layer.
#' 
#' @references
#' Cox, P.M., Betts, R.A., Bunton, C.B., Essery, R.L.H., Rowntree, P.R., and Smith, J. 1999. The impact of new land surface physics on the GCM simulation of climate and climate sensitivity. Climate Dynamics 15: 183–203.
#' 
#' Dharssi, I., Vidale, P.L., Verhoef, A., MacPherson, B., Jones, C., and Best, M. 2009. New soil physical properties implemented in the Unified Model at PS18. 9–12.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{soil}}
#' 
#' @examples
#' examplesoil = soil(defaultSoilParams())
#' soil_thermalConductivity(examplesoil)
#' soil_thermalCapacity(examplesoil)
#' 
#' #Values change when altering water content (drier layers have lower conductivity and capacity)
#' examplesoil$W = c(0.1, 0.4, 0.7, 1.0)
#' soil_thermalConductivity(examplesoil)
#' soil_thermalCapacity(examplesoil)
#' 
#' @name soil_thermodynamics
soil_thermalCapacity <- function(soil, model = "SX") {
    .Call(`_medfate_thermalCapacity`, soil, model)
}

#' @rdname soil_thermodynamics
soil_thermalConductivity <- function(soil, model = "SX") {
    .Call(`_medfate_thermalConductivity`, soil, model)
}

#' @name soil_thermodynamics
soil_temperatureGradient <- function(dVec, Temp) {
    .Call(`_medfate_temperatureGradient`, dVec, Temp)
}

#' @name soil_thermodynamics
soil_temperatureChange <- function(dVec, Temp, sand, clay, W, Theta_FC, Gdown) {
    .Call(`_medfate_temperatureChange`, dVec, Temp, sand, clay, W, Theta_FC, Gdown)
}

#' Single-day simulation
#'
#' Function \code{spwb_day} performs water balance for a single day and \code{growth_day} 
#' performs water and carbon balance for a single day.
#' 
#' @param x An object of class \code{\link{spwbInput}} or \code{\link{growthInput}}.
#' @param date Date as string "yyyy-mm-dd".
#' @param meteovec A named numerical vector with weather data. See variable names in parameter \code{meteo} of \code{\link{spwb}}.
#' @param latitude Latitude (in degrees).
#' @param elevation,slope,aspect Elevation above sea level (in m), slope (in degrees) and aspect (in degrees from North). 
#' @param runon Surface water amount running on the target area from upslope (in mm).
#' @param modifyInput Boolean flag to indicate that the input \code{x} object is allowed to be modified during the simulation.
#' 
#' @details
#' The simulation functions allow using three different sub-models of transpiration and photosynthesis:
#' \itemize{
#'   \item{The sub-model corresponding to 'Granier' transpiration mode is illustrated by function \code{\link{transp_transpirationGranier}} and was described in De Caceres et al. (2015),
#'   and implements an approach originally described in Granier et al. (1999).} 
#'   \item{The sub-model corresponding to 'Sperry' transpiration mode is illustrated by function \code{\link{transp_transpirationSperry}} and was described in De Caceres et al. (2021), and
#'   implements a modelling approach originally described in Sperry et al. (2017).}  
#'   \item{The sub-model corresponding to 'Cochard' transpiration mode is illustrated by function \code{\link{transp_transpirationCochard}} and was described for model SurEau-Ecos v2.0 in Ruffault et al. (2022).} 
#' }
#' Simulations using the 'Sperry' or 'Cochard' transpiration mode are computationally much more expensive than 'Granier'.
#' 
#' @return
#' Function \code{spwb_day()} returns a list of class \code{spwb_day} with the 
#' following elements:
#' \itemize{
#'   \item{\code{"cohorts"}: A data frame with cohort information, copied from \code{\link{spwbInput}}.}
#'   \item{\code{"topography"}: Vector with elevation, slope and aspect given as input.} 
#'   \item{\code{"weather"}: A vector with the input weather.}
#'   \item{\code{"WaterBalance"}: A vector of water balance components (rain, snow, net rain, infiltration, ...) for the simulated day, equivalent to one row of 'WaterBalance' object given in \code{\link{spwb}}.}
#'   \item{\code{"Soil"}: A data frame with results for each soil layer:
#'     \itemize{
#'       \item{\code{"SoilEvaporation"}: Water evaporated from the soil surface (in mm).}
#'       \item{\code{"HydraulicInput"}: Water entering each soil layer from other layers, transported via plant hydraulic network (in mm) (only for \code{transpirationMode = "Sperry"}).}
#'       \item{\code{"HydraulicOutput"}: Water leaving each soil layer (going to other layers or the transpiration stream) (in mm) (only for \code{transpirationMode = "Sperry"}).}
#'       \item{\code{"PlantExtraction"}: Water extracted by plants from each soil layer (in mm).}
#'       \item{\code{"psi"}: Soil water potential (in MPa).}
#'     }
#'   }
#'   \item{\code{"Stand"}: A named vector with with stand values for the simulated day, equivalent to one row of 'Stand' object returned by \code{\link{spwb}}.}
#'   \item{\code{"Plants"}: A data frame of results for each plant cohort (see \code{\link{transp_transpirationGranier}} or \code{\link{transp_transpirationSperry}}).}
#' }
#' The following items are only returned when \code{transpirationMode = "Sperry"} or  \code{transpirationMode = "Cochard"}:
#' \itemize{
#'   \item{\code{"EnergyBalance"}: Energy balance of the stand (see \code{\link{transp_transpirationSperry}}).}
#'   \item{\code{"RhizoPsi"}: Minimum water potential (in MPa) inside roots, after crossing rhizosphere, per cohort and soil layer.}
#'   \item{\code{"SunlitLeaves"} and \code{"ShadeLeaves"}: For each leaf type, a data frame with values of LAI, Vmax298 and Jmax298 for leaves of this type in each plant cohort.}
#'   \item{\code{"ExtractionInst"}: Water extracted by each plant cohort during each time step.}
#'   \item{\code{"PlantsInst"}: A list with instantaneous (per time step) results for each plant cohort (see \code{\link{transp_transpirationSperry}}).}
#'   \item{\code{"LightExtinction"}: A list of information regarding radiation balance through the canopy, as returned by function \code{\link{light_instantaneousLightExtinctionAbsortion}}.}
#'   \item{\code{"CanopyTurbulence"}: Canopy turbulence (see \code{\link{wind_canopyTurbulence}}).}
#' }
#'   
#' @references
#' De \enc{Cáceres}{Caceres} M, \enc{Martínez}{Martinez}-Vilalta J, Coll L, Llorens P, Casals P, Poyatos R, Pausas JG, Brotons L. (2015) Coupling a water balance model with forest inventory data to predict drought stress: the role of forest structural changes vs. climate changes. Agricultural and Forest Meteorology 213: 77-90 (doi:10.1016/j.agrformet.2015.06.012).
#' 
#' De \enc{Cáceres}{Caceres} M, Mencuccini M, Martin-StPaul N, Limousin JM, Coll L, Poyatos R, Cabon A, Granda V, Forner A, Valladares F, \enc{Martínez}{Martinez}-Vilalta J (2021) Unravelling the effect of species mixing on water use and drought stress in holm oak forests: a modelling approach. Agricultural and Forest Meteorology 296 (doi:10.1016/j.agrformet.2020.108233).
#' 
#' Granier A, \enc{Bréda}{Breda} N, Biron P, Villette S (1999) A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands. Ecol Modell 116:269–283. https://doi.org/10.1016/S0304-3800(98)00205-1.
#' 
#' Ruffault J, Pimont F, Cochard H, Dupuy JL, Martin-StPaul N (2022) 
#' SurEau-Ecos v2.0: a trait-based plant hydraulics model for simulations of plant water status and drought-induced mortality at the ecosystem level.
#' Geoscientific Model Development 15, 5593-5626 (doi:10.5194/gmd-15-5593-2022).
#' 
#' Sperry, J. S., M. D. Venturas, W. R. L. Anderegg, M. Mencuccini, D. S. Mackay, Y. Wang, and D. M. Love. 2017. Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost. Plant Cell and Environment 40, 816-830 (doi: 10.1111/pce.12852).
#' 
#' @author
#' \itemize{
#'   \item{Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF}
#'   \item{Nicolas Martin-StPaul, URFM-INRAE}
#' }
#' 
#' @seealso
#' \code{\link{spwbInput}}, \code{\link{spwb}},  \code{\link{plot.spwb_day}},  
#' \code{\link{growthInput}}, \code{\link{growth}},  \code{\link{plot.growth_day}}  
#' 
#' @examples
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#' 
#' # Day to be simulated
#' d <- 100
#' meteovec <- unlist(examplemeteo[d,])
#' date <- rownames(examplemeteo)[d]
#' 
#' #Simulate water balance one day only (Granier mode)
#' examplesoil <- soil(defaultSoilParams(4))
#' x1 <- forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' sd1 <- spwb_day(x1, date, meteovec,  
#'                 latitude = 41.82592, elevation = 100, slope=0, aspect=0) 
#' 
#' #Simulate water balance for one day only (Sperry mode)
#' control <- defaultControl("Sperry")
#' x2 <- forest2spwbInput(exampleforestMED, examplesoil, SpParamsMED, control)
#' sd2 <-spwb_day(x2, date, meteovec,
#'               latitude = 41.82592, elevation = 100, slope=0, aspect=0)
#' 
#' #Plot plant transpiration (see function 'plot.swb.day()')
#' plot(sd2)
#' 
#' #Simulate water balance for one day only (Cochard mode)
#' control <- defaultControl("Cochard")
#' x3 <- forest2spwbInput(exampleforestMED, examplesoil, SpParamsMED, control)
#' sd3 <-spwb_day(x3, date, meteovec,
#'               latitude = 41.82592, elevation = 100, slope=0, aspect=0)
#' 
#' 
#' #Simulate water and carbon balance for one day only (Granier mode)
#' control <- defaultControl("Granier")
#' x4  <- forest2growthInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' sd4 <- growth_day(x4, date, meteovec,
#'                 latitude = 41.82592, elevation = 100, slope=0, aspect=0)
#' 
#' @name spwb_day
spwb_day <- function(x, date, meteovec, latitude, elevation, slope, aspect, runon = 0.0, modifyInput = TRUE) {
    .Call(`_medfate_spwbDay`, x, date, meteovec, latitude, elevation, slope, aspect, runon, modifyInput)
}

#' Soil-plant water balance
#' 
#' Function \code{spwb()} is a water balance model that determines changes in soil moisture, 
#' soil water potentials, plant transpiration and drought stress at daily steps for a given forest stand 
#' during a period specified in the input climatic data. Function \code{pwb()} performs plant water balance 
#' only (i.e. soil moisture dynamics is an input) at daily steps for a given forest stand 
#' during a period specified in the input climatic data. On both simulation functions plant transpiration 
#' and photosynthesis processes are conducted with different level of detail depending on the transpiration mode.
#' 
#' @param x An object of class \code{\link{spwbInput}}.
#' @param meteo A data frame with daily meteorological data series. 
#' Row names of the data frame should correspond to date strings with format "yyyy-mm-dd" (see \code{\link{Date}}). Alternatively,
#' a column \code{dates} can contain \code{\link{Date}} or \code{\link{POSIXct}} classes.
#' The following columns are required and cannot have missing values:
#'   \itemize{
#'     \item{\code{MinTemperature}: Minimum temperature (in degrees Celsius).}
#'     \item{\code{MaxTemperature}: Maximum temperature (in degrees Celsius).}
#'     \item{\code{Precipitation}: Precipitation (in mm).}
#'   }
#' The following columns are required but can contain missing values (NOTE: missing values will raise warnings):
#'   \itemize{
#'     \item{\code{MinRelativeHumidity}: Minimum relative humidity (in percent).}
#'     \item{\code{MaxRelativeHumidity}: Maximum relative humidity (in percent).}
#'     \item{\code{Radiation}: Solar radiation (in MJ/m2/day).}
#'   }
#' The following columns are optional:
#'   \itemize{
#'     \item{\code{WindSpeed}: Above-canopy wind speed (in m/s). This column may not exist, or can be left with \code{NA} values. In both cases simulations will assume a constant value specified in \code{\link{defaultControl}}.}
#'     \item{\code{CO2}: Atmospheric (above-canopy) CO2 concentration (in ppm). This column may not exist, or can be left with \code{NA} values. In both cases simulations will assume a constant value specified in \code{\link{defaultControl}}.}
#'     \item{\code{Patm}: Atmospheric pressure (in kPa). This column may not exist, or can be left with \code{NA} values. In both cases, a value is estimated from elevation.}
#'   }
#' @param latitude Latitude (in degrees).
#' @param elevation,slope,aspect Elevation above sea level (in m), slope (in degrees) and aspect (in degrees from North).
#' @param CO2ByYear A named numeric vector with years as names and atmospheric CO2 concentration (in ppm) as values. Used to specify annual changes in CO2 concentration along the simulation (as an alternative to specifying daily values in \code{meteo}).
#' 
#' @details 
#' The simulation functions allow using three different sub-models of transpiration and photosynthesis:
#' \itemize{
#'   \item{The sub-model corresponding to 'Granier' transpiration mode is illustrated by function \code{\link{transp_transpirationGranier}} and was described in De Caceres et al. (2015),
#'   and implements an approach originally described in Granier et al. (1999).} 
#'   \item{The sub-model corresponding to 'Sperry' transpiration mode is illustrated by function \code{\link{transp_transpirationSperry}} and was described in De Caceres et al. (2021), and
#'   implements a modelling approach originally described in Sperry et al. (2017).}  
#'   \item{The sub-model corresponding to 'Cochard' transpiration mode is illustrated by function \code{\link{transp_transpirationCochard}} and was described for model SurEau-Ecos v2.0 in Ruffault et al. (2022).} 
#' }
#' Simulations using the 'Sperry' or 'Cochard' transpiration mode are computationally much more expensive than 'Granier'.
#' 
#' @return
#' Function \code{spwb} returns a list of class 'spwb' whereas function \code{pwb} returns a list of class 'pwb'. 
#' There are many elements in common in these lists, so they are listed here together:
#' \itemize{
#'   \item{\code{"latitude"}: Latitude (in degrees) given as input.} 
#'   \item{\code{"topography"}: Vector with elevation, slope and aspect given as input.} 
#'   \item{\code{"weather"}: A copy of the input weather data frame.}
#'   \item{\code{"spwbInput"}: An copy of the object \code{x} of class \code{\link{spwbInput}} given as input.}
#'   \item{\code{"spwbOutput"}: An copy of the final state of the object \code{x} of class \code{\link{spwbInput}}.}
#'   \item{\code{"WaterBalance"}: A data frame where different variables (in columns) are given for each simulated day (in rows):}
#'   \itemize{
#'     \item{\code{"PET"}: Potential evapotranspiration (in mm).}
#'     \item{\code{"Precipitation"}: Input precipitation (in mm).}
#'     \item{\code{"Rain"}: Precipitation as rain (in mm).}
#'     \item{\code{"Snow"}: Precipitation as snow (in mm).}
#'     \item{\code{"NetRain"}: Net rain, after accounting for interception (in mm).}
#'     \item{\code{"Infiltration"}: The amount of water infiltrating into the soil (in mm).}
#'     \item{\code{"Runoff"}: The amount of water exported via surface runoff (in mm).}
#'     \item{\code{"DeepDrainage"}: The amount of water exported via deep drainage (in mm).}
#'     \item{\code{"Evapotranspiration"}: Evapotranspiration (in mm).}
#'     \item{\code{"SoilEvaporation"}: Bare soil evaporation (in mm).}
#'     \item{\code{"HerbTranspiration"}: Transpiration due to the herbaceous layer (in mm).}
#'     \item{\code{"PlantExtraction"}: Amount of water extracted from soil by woody plants (in mm).}
#'     \item{\code{"Transpiration"}: Woody plant transpiration (in mm).}
#'     \item{\code{"HydraulicRedistribution"}: Water redistributed among soil layers, transported through the plant hydraulic network.}
#'   }
#'   \item{\code{"EnergyBalance"}: A data frame with the daily values of energy balance components for the soil and the canopy (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}).}
#'   \item{\code{"Temperature"}: A data frame with the daily values of minimum/mean/maximum temperatures for the atmosphere (input), canopy and soil (only for \code{transpirationMode = "Sperry"} or \code{transpirationMode = "Cochard"}).}
#'   \item{\code{"Soil"}: A data frame where different variables (in columns) are given for each simulated day (in rows):}
#'   \itemize{
#'     \item{\code{"W.1"}, \code{...}, \code{"W.k"}: Relative soil moisture content (relative to field capacity) in each soil layer.}
#'     \item{\code{"ML.1"}, \code{...}, \code{"ML.k"}: Soil water volume in each soil layer (in L/m2).}
#'     \item{\code{"MLTot"}: Total soil water volume (in L/m2).}
#'     \item{\code{"SWE"}: Snow water equivalent (mm) of the snow pack.}
#'     \item{\code{"PlantExt.1"}, \code{...}, \code{"PlantExt.k"}: Plant extraction from each soil layer (in mm).}
#'     \item{\code{"HydraulicInput.1"}, \code{...}, \code{"HydraulicInput.k"}: Water that entered the layer coming from other layers and transported via the plant hydraulic network (in mm).}
#'     \item{\code{"psi.1"}, \code{...}, \code{"psi.k"}: Soil water potential in each soil layer (in MPa).}
#'   }
#'   \item{\code{"Stand"}: A data frame where different variables (in columns) are given for each simulated day (in rows):}
#'   \itemize{
#'     \item{\code{"LAI"}: LAI of the stand (including the herbaceous layer and live + dead leaves of woody plants) (in m2/m2).}
#'     \item{\code{"LAIherb"}: LAI of the herbaceous layer (in m2/m2).}
#'     \item{\code{"LAIlive"}: LAI of the woody plants assuming all leaves are unfolded (in m2/m2).}
#'     \item{\code{"LAIexpanded"}: LAI of the woody plants with leaves actually unfolded (in m2/m2).}
#'     \item{\code{"LAIdead"}: LAI of the woody plants corresponding to dead leaves (in m2/m2).}
#'     \item{\code{"Cm"}: Water retention capacity of the canopy (in mm) (accounting for leaf phenology).}
#'     \item{\code{"LgroundPAR"}: The percentage of PAR that reaches the ground (accounting for leaf phenology).}
#'     \item{\code{"LgroundSWR"}: The percentage of SWR that reaches the ground (accounting for leaf phenology).}
#'   }
#'   \item{\code{"Plants"}: A list of daily results for plant cohorts (see below).}
#'   \item{\code{"subdaily"}: A list of objects of class \code{\link{spwb_day}}, one per day simulated (only if required in \code{control} parameters, see \code{\link{defaultControl}}).}
#' }
#' 
#' When \code{transpirationMode = "Granier"}, element \code{"Plants"} is a list with the following subelements:
#'   \itemize{
#'     \item{\code{"LAI"}: A data frame with the daily leaf area index for each plant cohort.}
#'     \item{\code{"LAIlive"}: A data frame with the daily leaf area index for each plant cohort, assuming all leaves are unfolded (in m2/m2).}
#'     \item{\code{"FPAR"}: A data frame with the fraction of PAR at the canopy level of each plant cohort. }
#'     \item{\code{"AbsorbedSWRFraction"}: A data frame with the fraction of SWR absorbed by each plant cohort. }
#'     \item{\code{"Transpiration"}: A data frame with the amount of daily transpiration (in mm) for each plant cohort.}
#'     \item{\code{"GrossPhotosynthesis"}: A data frame with the amount of daily gross photosynthesis (in g C·m-2) for each plant cohort. }
#'     \item{\code{"PlantPsi"}: A data frame with the average daily water potential of each plant (in MPa).}
#'     \item{\code{"StemPLC"}: A data frame with the average daily proportion of stem conductance loss of each plant ([0-1]).}
#'     \item{\code{"PlantWaterBalance"}: A data frame with the daily balance between transpiration and soil water extraction for each plant cohort. }
#'     \item{\code{"LeafRWC"}: A data frame with the average daily leaf relative water content of each plant (in percent).}
#'     \item{\code{"StemRWC"}: A data frame with the average daily stem relative water content of each plant (in percent). }
#'     \item{\code{"LFMC"}: A data frame with the daily live fuel moisture content (in percent of dry weight).}
#'     \item{\code{"PlantStress"}: A data frame with the amount of daily stress [0-1] suffered by each plant cohort (relative whole-plant conductance).}
#'   }
#' If \code{transpirationMode="Sperry"} or \code{transpirationMode="Cochard"}, element \code{"Plants"} is a list with the following subelements:
#'   \itemize{
#'     \item{\code{"LAI"}: A data frame with the daily leaf area index for each plant cohort.}
#'     \item{\code{"AbsorbedSWR"}: A data frame with the daily SWR absorbed by each plant cohort.}
#'     \item{\code{"NetLWR"}: A data frame with the daily net LWR by each plant cohort.}
#'     \item{\code{"Transpiration"}: A data frame with the amount of daily transpiration (in mm) for each plant cohorts.}
#'     \item{\code{"GrossPhotosynthesis"}: A data frame with the amount of daily gross photosynthesis (in g C·m-2) for each plant cohort. }
#'     \item{\code{"NetPhotosynthesis"}: A data frame with the amount of daily net photosynthesis (in g C·m-2) for each plant cohort. }
#'     \item{\code{"dEdP"}: A data frame with mean daily values of soil-plant conductance (derivative of the supply function) for each plant cohort.}
#'     \item{\code{"PlantWaterBalance"}: A data frame with the daily balance between transpiration and soil water extraction for each plant cohort. }
#'     \item{\code{"SunlitLeaves"} and \code{"ShadeLeaves"}: A list with daily results for sunlit and shade leaves:
#'       \itemize{
#'         \item{\code{"PsiMin"}: A data frame with the minimum (midday) daily sunlit or shade leaf water potential (in MPa). }
#'         \item{\code{"PsiMax"}: A data frame with the maximum (predawn) daily sunlit or shade leaf water potential (in MPa). }
#'       }
#'     }
#'     \item{\code{"LeafPsiMin"}: A data frame with the minimum (midday) daily (average) leaf water potential of each plant (in MPa).}
#'     \item{\code{"LeafPsiMax"}: A data frame with the maximum (predawn) daily (average) leaf water potential of each plant (in MPa).}
#'     \item{\code{"LeafRWC"}: A data frame with the average daily leaf relative water content of each plant (in percent).}
#'     \item{\code{"StemRWC"}: A data frame with the average daily stem relative water content of each plant (in percent). }
#'     \item{\code{"LFMC"}: A data frame with the daily live fuel moisture content (in percent of dry weight).}
#'     \item{\code{"StemPsi"}: A data frame with the minimum daily stem water potential of each plant (in MPa). }
#'     \item{\code{"StemPLC"}: A data frame with the average daily proportion of stem conductance loss of each plant ([0-1]).}
#'     \item{\code{"RootPsi"}: A data frame with the minimum daily root water potential of each plant (in MPa). }
#'     \item{\code{"RhizoPsi"}: A list of data frames (one per plant cohort) with the minimum daily root water potential of each plant (in MPa).}
#'     \item{\code{"PlantStress"}: A data frame with the amount of daily stress [0-1] suffered by each plant cohort (relative whole-plant conductance).}
#'   }
#' 
#' @references
#' De \enc{Cáceres}{Caceres} M, \enc{Martínez}{Martinez}-Vilalta J, Coll L, Llorens P, Casals P, Poyatos R, Pausas JG, Brotons L. (2015) Coupling a water balance model with forest inventory data to predict drought stress: the role of forest structural changes vs. climate changes. Agricultural and Forest Meteorology 213: 77-90 (doi:10.1016/j.agrformet.2015.06.012).
#' 
#' De \enc{Cáceres}{Caceres} M, Mencuccini M, Martin-StPaul N, Limousin JM, Coll L, Poyatos R, Cabon A, Granda V, Forner A, Valladares F, \enc{Martínez}{Martinez}-Vilalta J (2021) Unravelling the effect of species mixing on water use and drought stress in holm oak forests: a modelling approach. Agricultural and Forest Meteorology 296 (doi:10.1016/j.agrformet.2020.108233).
#' 
#' Granier A, \enc{Bréda}{Breda} N, Biron P, Villette S (1999) A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands. Ecol Modell 116:269–283. https://doi.org/10.1016/S0304-3800(98)00205-1.
#' 
#' Ruffault J, Pimont F, Cochard H, Dupuy JL, Martin-StPaul N (2022) 
#' SurEau-Ecos v2.0: a trait-based plant hydraulics model for simulations of plant water status and drought-induced mortality at the ecosystem level.
#' Geoscientific Model Development 15, 5593-5626 (doi:10.5194/gmd-15-5593-2022).
#' 
#' Sperry, J. S., M. D. Venturas, W. R. L. Anderegg, M. Mencuccini, D. S. Mackay, Y. Wang, and D. M. Love. 2017. Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost. Plant Cell and Environment 40, 816-830 (doi: 10.1111/pce.12852).
#' 
#' @author
#' \itemize{
#'   \item{Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF}
#'   \item{Nicolas Martin-StPaul, URFM-INRAE}
#' }
#' 
#' @seealso 
#' \code{\link{spwbInput}}, \code{\link{spwb_day}}, \code{\link{plot.spwb}}, 
#' \code{\link{forest}}
#' 
#' @examples
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Initialize soil with default soil params (4 layers)
#' examplesoil <- soil(defaultSoilParams(4))
#' 
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#' 
#' #Initialize input
#' x1 <- forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' S1 <- spwb(x1, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' #Plot results
#' plot(S1)
#' 
#' #Monthly summary (averages) of soil water balance
#' summary(S1, freq="months",FUN=mean, output="Soil")
#'                   
#' \donttest{
#' #Switch to 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' 
#' #Initialize input
#' x2 <- forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' S2 <- spwb(x2, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' #Switch to 'Cochard' transpiration mode
#' control <- defaultControl("Cochard")
#' 
#' #Initialize input
#' x3 <- forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' S3 <- spwb(x3, examplemeteo, latitude = 41.82592, elevation = 100)
#' }
#'                 
#' @name spwb
spwb <- function(x, meteo, latitude, elevation = NA_real_, slope = NA_real_, aspect = NA_real_, CO2ByYear = numeric(0)) {
    .Call(`_medfate_spwb`, x, meteo, latitude, elevation, slope, aspect, CO2ByYear)
}

#' @rdname spwb
#' 
#' @param W A matrix with the same number of rows as \code{meteo} and as many columns as soil layers, containing the soil moisture of each layer as proportion of field capacity.
#' @param canopyEvaporation A vector of daily canopy evaporation (from interception) values (mm). The length should match the number of rows in \code{meteo}.
#' @param snowMelt A vector of daily snow melt values (mm). The length should match the number of rows in \code{meteo}.
#' @param soilEvaporation A vector of daily bare soil evaporation values (mm). The length should match the number of rows in \code{meteo}.
#' @param herbTranspiration A vector of daily herbaceous transpiration values (mm). The length should match the number of rows in \code{meteo}.
#' 
pwb <- function(x, meteo, W, latitude, elevation = NA_real_, slope = NA_real_, aspect = NA_real_, canopyEvaporation = numeric(0), snowMelt = numeric(0), soilEvaporation = numeric(0), herbTranspiration = numeric(0), CO2ByYear = numeric(0)) {
    .Call(`_medfate_pwb`, x, meteo, W, latitude, elevation, slope, aspect, canopyEvaporation, snowMelt, soilEvaporation, herbTranspiration, CO2ByYear)
}

#' Tissue moisture functions
#' 
#' Set of functions used to calculate tissue moisture from water potential and viceversa.
#' 
#' @param psiSym,psiApo Symplastic or apoplastic water potential (MPa).
#' @param RWC Relative water content [0-1].
#' @param pi0 Full turgor osmotic potential (MPa).
#' @param epsilon Bulk modulus of elasticity (MPa).
#' @param c,d Parameters of the xylem vulnerability curve.
#' @param af Apoplastic fraction (proportion) in the segment (e.g. leaf or stem).
#' @param L Vector with the length of coarse roots (mm) for each soil layer.
#' @param V Vector with the proportion [0-1] of fine roots within each soil layer.
#' @param Al2As Leaf area to sapwood area (in m2·m-2).
#' @param height Plant height (in cm).
#' @param SLA Specific leaf area (mm2·mg-1).
#' @param wd Wood density (g·cm-3).
#' @param ld Leaf tissue density (g·cm-3).
#' 
#' @return
#' Values returned for each function are:
#' \itemize{
#'   \item{\code{moisture_symplasticRWC}: Relative water content [0-1] of the symplastic fraction.}
#'   \item{\code{moisture_apoplasticRWC}: Relative water content [0-1] of the apoplastic fraction.}
#'   \item{\code{moisture_symplasticWaterPotential}: Water potential (in MPa) of the symplastic fraction.}
#'   \item{\code{moisture_apoplasticWaterPotential}: Water potential (in MPa) of the apoplastic fraction.}
#'   \item{\code{moisture_turgorLossPoint}: Water potential (in MPa) corresponding to turgor loss point.}
#'   \item{\code{moisture_segmentRWC}: Segment relative water content [0-1].}
#'   \item{\code{water_plant}: A vector of water content (mm) per plant cohort.}
#' }
#' 
#' @references
#' Bartlett, M.K., Scoffoni, C., Sack, L. 2012. The determinants of leaf turgor loss point and prediction of drought tolerance of species and biomes: a global meta-analysis. Ecology Letters 15: 393–405.
#' 
#' \enc{Hölttä}{Holtta}, T., Cochard, H., Nikinmaa, E., Mencuccini, M. 2009. Capacitive effect of cavitation in xylem conduits: Results from a dynamic model. Plant, Cell and Environment 32: 10–21.
#' 
#' Martin-StPaul, N., Delzon, S., Cochard, H. 2017. Plant resistance to drought depends on timely stomatal closure. Ecology Letters 20: 1437–1447.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{hydraulics_psi2K}}, \code{\link{hydraulics_supplyFunctionPlot}}, \code{\link{spwb}}, \code{\link{soil}}
#' 
#' @examples
#' psi = seq(-10,0, by=0.1)
#' rwc_s = rep(NA, length(psi))
#' for(i in 1:length(psi)) rwc_s[i] = moisture_symplasticRWC(psi[i],-3,12)
#' plot(psi, rwc_s, type="l", xlab="Water potential (MPa)", ylab = "Symplasmic RWC")
#' 
#' @name moisture
moisture_sapwoodWaterCapacity <- function(Al2As, height, V, L, wd) {
    .Call(`_medfate_sapwoodWaterCapacity`, Al2As, height, V, L, wd)
}

#' @rdname moisture
moisture_leafWaterCapacity <- function(SLA, ld) {
    .Call(`_medfate_leafWaterCapacity`, SLA, ld)
}

#' @rdname moisture
moisture_turgorLossPoint <- function(pi0, epsilon) {
    .Call(`_medfate_turgorLossPoint`, pi0, epsilon)
}

#' @rdname moisture
moisture_symplasticRWC <- function(psiSym, pi0, epsilon) {
    .Call(`_medfate_symplasticRelativeWaterContent`, psiSym, pi0, epsilon)
}

#' @rdname moisture
moisture_symplasticPsi <- function(RWC, pi0, epsilon) {
    .Call(`_medfate_symplasticWaterPotential`, RWC, pi0, epsilon)
}

#' @rdname moisture
moisture_apoplasticRWC <- function(psiApo, c, d) {
    .Call(`_medfate_apoplasticRelativeWaterContent`, psiApo, c, d)
}

#' @rdname moisture
moisture_apoplasticPsi <- function(RWC, c, d) {
    .Call(`_medfate_apoplasticWaterPotential`, RWC, c, d)
}

#' @rdname moisture
moisture_tissueRWC <- function(psiSym, pi0, epsilon, psiApo, c, d, af) {
    .Call(`_medfate_tissueRelativeWaterContent`, psiSym, pi0, epsilon, psiApo, c, d, af)
}

#' @rdname moisture
plant_water <- function(x) {
    .Call(`_medfate_plantWaterContent`, x)
}

#' @rdname transp_modes
#' 
#' @param canopyEvaporation Canopy evaporation (from interception) for \code{day} (mm).
#' @param snowMelt Snow melt values  for \code{day} (mm).
#' @param soilEvaporation Bare soil evaporation for \code{day} (mm).
#' @param herbTranspiration Transpiration of herbaceous plants for \code{day} (mm).
#' @param stepFunctions An integer to indicate a simulation step for which photosynthesis and profit maximization functions are desired.
#' 
#' 
transp_transpirationSperry <- function(x, meteo, day, latitude, elevation, slope, aspect, canopyEvaporation = 0.0, snowMelt = 0.0, soilEvaporation = 0.0, herbTranspiration = 0.0, stepFunctions = NA_integer_, modifyInput = TRUE) {
    .Call(`_medfate_transpirationSperry`, x, meteo, day, latitude, elevation, slope, aspect, canopyEvaporation, snowMelt, soilEvaporation, herbTranspiration, stepFunctions, modifyInput)
}

#' @rdname transp_modes
transp_transpirationCochard <- function(x, meteo, day, latitude, elevation, slope, aspect, canopyEvaporation = 0.0, snowMelt = 0.0, soilEvaporation = 0.0, herbTranspiration = 0.0, modifyInput = TRUE) {
    .Call(`_medfate_transpirationCochard`, x, meteo, day, latitude, elevation, slope, aspect, canopyEvaporation, snowMelt, soilEvaporation, herbTranspiration, modifyInput)
}

#' Transpiration modes
#' 
#' High-level sub-models representing transpiration, plant hydraulics, photosynthesis and water relations 
#' within plants. 
#' 
#' Three sub-models are available: 
#' \itemize{
#'   \item{Sub-model in function \code{transp_transpirationGranier} was described in De \enc{Cáceres}{Caceres} et al. (2015), 
#'   and implements an approach originally described in Granier et al. (1999).} 
#'   \item{Sub-model in function \code{transp_transpirationSperry} was described in De \enc{Cáceres}{Caceres} et al. (2021), and
#'   implements a modelling approach originally described in Sperry et al. (2017).} 
#'   \item{Sub-model in function \code{transp_transpirationCochard} was described for SurEau-Ecos v2.0 model in Ruffault et al. (2022).} 
#' }
#' 
#' @param x An object of class \code{\link{spwbInput}} or \code{\link{growthInput}}, built using the 'Granier', 'Sperry' or 'Cochard' transpiration modes.
#' @param meteo A data frame with daily meteorological data series (see \code{\link{spwb}}).
#' @param day An integer to identify a day (row) within the \code{meteo} data frame.
#' @param latitude Latitude (in degrees).
#' @param elevation,slope,aspect Elevation above sea level (in m), slope (in degrees) and aspect (in degrees from North).
#' @param modifyInput Boolean flag to indicate that the input \code{x} object is allowed to be modified during the simulation.
#' 
#' @return
#' A list with the following elements:
#' \itemize{
#'   \item{\code{"cohorts"}: A data frame with cohort information, copied from \code{\link{spwbInput}}.}
#'   \item{\code{"Stand"}: A vector of stand-level variables.}
#'   \item{\code{"Plants"}: A data frame of results for each plant cohort. When using \code{transp_transpirationGranier}, element \code{"Plants"} includes:
#'     \itemize{
#'       \item{\code{"LAI"}: Leaf area index of the plant cohort.}
#'       \item{\code{"LAIlive"}: Leaf area index of the plant cohort, assuming all leaves are unfolded.}
#'       \item{\code{"AbsorbedSWRFraction"}: Fraction of SWR absorbed by each cohort.}
#'       \item{\code{"Transpiration"}: Transpirated water (in mm) corresponding to each cohort.}
#'       \item{\code{"GrossPhotosynthesis"}: Gross photosynthesis (in gC/m2) corresponding to each cohort.}
#'       \item{\code{"psi"}: Water potential (in MPa) of the plant cohort (average over soil layers).}
#'       \item{\code{"DDS"}: Daily drought stress [0-1] (relative whole-plant conductance).}
#'     }
#'   When using \code{transp_transpirationSperry} or \code{transp_transpirationCochard}, element \code{"Plants"} includes:
#'     \itemize{
#'       \item{\code{"LAI"}: Leaf area index of the plant cohort.}
#'       \item{\code{"LAIlive"}: Leaf area index of the plant cohort, assuming all leaves are unfolded.}
#'       \item{\code{"Extraction"}: Water extracted from the soil (in mm) for each cohort.}
#'       \item{\code{"Transpiration"}: Transpirated water (in mm) corresponding to each cohort.}
#'       \item{\code{"GrossPhotosynthesis"}: Gross photosynthesis (in gC/m2) corresponding to each cohort.}
#'       \item{\code{"NetPhotosynthesis"}: Net photosynthesis (in gC/m2) corresponding to each cohort.}
#'       \item{\code{"RootPsi"}: Minimum water potential (in MPa) at the root collar.}
#'       \item{\code{"StemPsi"}: Minimum water potential (in MPa) at the stem.}
#'       \item{\code{"StemPLC"}: Proportion of conductance loss in stem.}
#'       \item{\code{"LeafPsiMin"}: Minimum (predawn) water potential (in MPa) at the leaf (representing an average leaf).}
#'       \item{\code{"LeafPsiMax"}: Maximum (midday) water potential (in MPa) at the leaf (representing an average leaf).}
#'       \item{\code{"LeafPsiMin_SL"}: Minimum (predawn) water potential (in MPa) at sunlit leaves.}
#'       \item{\code{"LeafPsiMax_SL"}: Maximum (midday) water potential (in MPa) at sunlit leaves.}
#'       \item{\code{"LeafPsiMin_SH"}: Minimum (predawn) water potential (in MPa) at shade leaves.}
#'       \item{\code{"LeafPsiMax_SH"}: Maximum (midday) water potential (in MPa) at shade leaves.}
#'       \item{\code{"dEdP"}: Overall soil-plant conductance (derivative of the supply function).}
#'       \item{\code{"DDS"}: Daily drought stress [0-1] (relative whole-plant conductance).}
#'       \item{\code{"StemRWC"}: Relative water content of stem tissue (including symplasm and apoplasm).}
#'       \item{\code{"LeafRWC"}: Relative water content of leaf tissue (including symplasm and apoplasm).}
#'       \item{\code{"LFMC"}: Live fuel moisture content (in percent of dry weight).}
#'       \item{\code{"WaterBalance"}: Plant water balance (extraction - transpiration).}
#'     }
#'   }
#'   \item{\code{"Extraction"}: A data frame with mm of water extracted from each soil layer (in columns) by each cohort (in rows).}
#' 
#'   The remaining items are only given by \code{transp_transpirationSperry} or \code{transp_transpirationCochard}:
#'   \item{\code{"EnergyBalance"}: A list with the following elements:
#'     \itemize{
#'       \item{\code{"Temperature"}: A data frame with the temperature of the atmosphere ('Tatm'), canopy ('Tcan') and soil ('Tsoil.1', 'Tsoil.2', ...) for each time step.}
#'       \item{\code{"CanopyEnergyBalance"}: A data frame with the components of the canopy energy balance (in W/m2) for each time step.}
#'       \item{\code{"SoilEnergyBalance"}: A data frame with the components of the soil energy balance (in W/m2) for each time step.}
#'     }  
#'   }
#'   \item{\code{"RhizoPsi"}: Minimum water potential (in MPa) inside roots, after crossing rhizosphere, per cohort and soil layer.}
#'   \item{\code{"Sunlitleaves"} and \code{"ShadeLeaves"}: Data frames for sunlit leaves and shade leaves and the following columns per cohort:
#'     \itemize{
#'       \item{\code{"LAI"}: Cumulative leaf area index of sunlit/shade leaves.}
#'       \item{\code{"Vmax298"}: Average maximum carboxilation rate for sunlit/shade leaves.}
#'       \item{\code{"Jmax298"}: Average maximum electron transport rate for sunlit/shade leaves.}
#'     }  
#'   }
#'   \item{\code{"ExtractionInst"}: Water extracted by each plant cohort during each time step.}
#'   \item{\code{"PlantsInst"}: A list with instantaneous (per time step) results for each plant cohort:
#'     \itemize{
#'       \item{\code{"E"}: A data frame with the cumulative transpiration (mm) for each plant cohort during each time step. }
#'       \item{\code{"Ag"}: A data frame with the cumulative gross photosynthesis (gC/m2) for each plant cohort during each time step. }
#'       \item{\code{"An"}: A data frame with the cumulative net photosynthesis (gC/m2) for each plant cohort during each time step. }
#'       \item{\code{"Sunlitleaves"} and \code{"ShadeLeaves"}: Lists with instantaneous (for each time step) results for sunlit leaves and shade leaves and the following items:
#'         \itemize{
#'           \item{\code{"Abs_SWR"}: A data frame with instantaneous absorbed short-wave radiation (SWR).} 
#'           \item{\code{"Net_LWR"}: A data frame with instantaneous net long-wave radiation (LWR).} 
#'           \item{\code{"An"}: A data frame with instantaneous net photosynthesis (in micromol/m2/s). }
#'           \item{\code{"Ci"}: A data frame with instantaneous intercellular CO2 concentration (in ppm). }
#'           \item{\code{"GW"}: A data frame with instantaneous stomatal conductance (in mol/m2/s). }
#'           \item{\code{"VPD"}: A data frame with instantaneous vapour pressure deficit (in kPa). }
#'           \item{\code{"Temp"}: A data frame with leaf temperature (in degrees Celsius). }
#'           \item{\code{"Psi"}: A data frame with leaf water potential (in MPa). }
#'         }
#'       }
#'       \item{\code{"dEdP"}: A data frame with the slope of the plant supply function (an estimation of whole-plant conductance).}
#'       \item{\code{"RootPsi"}: A data frame with root crown water potential (in MPa) for each plant cohort during each time step.}
#'       \item{\code{"StemPsi"}: A data frame with stem water potential (in MPa) for each plant cohort during each time step.}
#'       \item{\code{"LeafPsi"}: A data frame with leaf (average) water potential (in MPa) for each plant cohort during each time step. }
#'       \item{\code{"StemPLC"}: A data frame with the proportion loss of conductance [0-1] for each plant cohort during each time step. }
#'       \item{\code{"StemRWC"}: A data frame with the (average) relative water content of stem tissue [0-1] for each plant cohort during each time step. }
#'       \item{\code{"LeafRWC"}: A data frame with the relative water content of leaf tissue [0-1] for each plant cohort during each time step. }
#'       \item{\code{"StemSympRWC"}: A data frame with the (average) relative water content of symplastic stem tissue [0-1] for each plant cohort during each time step. }
#'       \item{\code{"LeafSympRWC"}: A data frame with the relative water content of symplastic leaf tissue [0-1] for each plant cohort during each time step. }
#'       \item{\code{"PWB"}: A data frame with plant water balance (extraction - transpiration).}
#'     }
#'   }
#'   \item{\code{"LightExtinction"}: A list of information regarding radiation balance through the canopy, as returned by function \code{\link{light_instantaneousLightExtinctionAbsortion}}.}
#'   \item{\code{"CanopyTurbulence"}: Canopy turbulence (see \code{\link{wind_canopyTurbulence}}).}
#'   \item{\code{"SupplyFunctions"}: If \code{stepFunctions} is not missing, a list of supply functions, photosynthesis functions and profit maximization functions.}
#' }
#' 
#' @references
#' De \enc{Cáceres}{Caceres} M, \enc{Martínez}{Martinez}-Vilalta J, Coll L, Llorens P, Casals P, Poyatos R, Pausas JG, Brotons L. (2015) Coupling a water balance model with forest inventory data to predict drought stress: the role of forest structural changes vs. climate changes. Agricultural and Forest Meteorology 213: 77-90 (doi:10.1016/j.agrformet.2015.06.012).
#' 
#' De \enc{Cáceres}{Caceres} M, Mencuccini M, Martin-StPaul N, Limousin JM, Coll L, Poyatos R, Cabon A, Granda V, Forner A, Valladares F, \enc{Martínez}{Martinez}-Vilalta J (2021) Unravelling the effect of species mixing on water use and drought stress in holm oak forests: a modelling approach. Agricultural and Forest Meteorology 296 (doi:10.1016/j.agrformet.2020.108233).
#' 
#' Granier A, \enc{Bréda}{Breda} N, Biron P, Villette S (1999) A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands. Ecol Modell 116:269–283. https://doi.org/10.1016/S0304-3800(98)00205-1.
#' 
#' Ruffault J, Pimont F, Cochard H, Dupuy JL, Martin-StPaul N (2022) 
#' SurEau-Ecos v2.0: a trait-based plant hydraulics model for simulations of plant water status and drought-induced mortality at the ecosystem level.
#' Geoscientific Model Development 15, 5593-5626 (doi:10.5194/gmd-15-5593-2022).
#' 
#' Sperry, J. S., M. D. Venturas, W. R. L. Anderegg, M. Mencuccini, D. S. Mackay, Y. Wang, and D. M. Love. 2017. Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost. Plant Cell and Environment 40, 816-830 (doi: 10.1111/pce.12852).
#' 
#' @author 
#' \itemize{
#'   \item{Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF}
#'   \item{Nicolas Martin-StPaul, URFM-INRAE}
#' }
#' 
#' @seealso \code{\link{spwb_day}}, \code{\link{plot.spwb_day}}
#' 
#' @examples
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Initialize soil with default soil params (4 layers)
#' examplesoil <- soil(defaultSoilParams(4))
#' 
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#' 
#' #Initialize input
#' x1 <- forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' # Transpiration according to Granier's model, plant water potential 
#' # and plant stress for a given day
#' t1 <- transp_transpirationGranier(x1, examplemeteo, 1, 
#'                                  latitude = 41.82592, elevation = 100, slope = 0, aspect = 0, 
#'                                  modifyInput = FALSE)
#' 
#' #Switch to 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' 
#' #Initialize input
#' x2 <- forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' # Transpiration according to Sperry's model
#' t2 <- transp_transpirationSperry(x2, examplemeteo, 1, 
#'                                 latitude = 41.82592, elevation = 100, slope = 0, aspect = 0,
#'                                 modifyInput = FALSE)
#'                                 
#' #Switch to 'Cochard' transpiration mode
#' control <- defaultControl("Cochard")
#' 
#' #Initialize input
#' x3 <- forest2spwbInput(exampleforestMED,examplesoil, SpParamsMED, control)
#' 
#' # Transpiration according to Cochard's model
#' t3 <- transp_transpirationCochard(x3, examplemeteo, 1, 
#'                                   latitude = 41.82592, elevation = 100, slope = 0, aspect = 0,
#'                                   modifyInput = FALSE)
#'                                 
#' @name transp_modes
transp_transpirationGranier <- function(x, meteo, day, latitude, elevation, slope, aspect, modifyInput = TRUE) {
    .Call(`_medfate_transpirationGranier`, x, meteo, day, latitude, elevation, slope, aspect, modifyInput)
}

#' Models for canopy turbulence
#' 
#' Models for canopy turbulence by Katul et al (2004).
#' 
#' @param zm A numeric vector with height values (m).
#' @param Cx Effective drag = Cd x leaf area density.
#' @param hm Canopy height (m).
#' @param d0 Zero displacement height (m).
#' @param z0 Momentum roughness height (m).
#' @param zmid A numeric vector of mid-point heights (in cm) for canopy layers.
#' @param LAD A numeric vector of leaf area density values (m3/m2).
#' @param canopyHeight Canopy height (in cm).
#' @param u Measured wind speed (m/s).
#' @param windMeasurementHeight Height of wind speed measurement with respect to canopy height (cm).
#' @param model Closure model.
#' 
#' @return 
#' Function \code{wind_canopyTurbulenceModel} returns a data frame of vertical profiles for variables:
#' \itemize{
#'   \item{\code{z1}: Height values.}
#'   \item{\code{U1}: U/u*, where U is mean velocity and u* is friction velocity.}
#'   \item{\code{dU1}: dUdz/u*, where dUdz is mean velocity gradient and u* is friction velocity.}
#'   \item{\code{epsilon1}: epsilon/(u^3/h) where epsilon is the turbulent kinetic dissipation rate, u* is friction velocity and h is canopy height.}
#'   \item{\code{k1}: k/(u*^2), where k is the turbulent kinetic energy and u* is friction velocity.}
#'   \item{\code{uw1}: <uw>/(u*^2), where <uw> is the Reynolds stress and u* is friction velocity.}
#'   \item{\code{Lmix1}: Mixing length.}
#' }
#' 
#' Function \code{wind_canopyTurbulence} returns a data frame of vertical profiles for transformed variables:
#'   \itemize{
#'     \item{\code{zmid}: Input mid-point heights (in cm) for canopy layers.}
#'     \item{\code{u}: Wind speed (m/s).}
#'     \item{\code{du}: Mean velocity gradient (1/s).}
#'     \item{\code{epsilon}: Turbulent kinetic dissipation rate.}
#'     \item{\code{k}: Turbulent kinetic energy.}
#'     \item{\code{uw}: Reynolds stress.}
#'   }
#' 
#' @details
#' Implementation in Rcpp of the K-epsilon canopy turbulence models by Katul et al (2004) originally in Matlab code (https://nicholas.duke.edu/people/faculty/katul/k_epsilon_model.htm).
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @references 
#' Katul GG, Mahrt L, Poggi D, Sanz C (2004) One- and two-equation models for canopy turbulence. Boundary-Layer Meteorol 113:81–109. https://doi.org/10.1023/B:BOUN.0000037333.48760.e5
#' 
#' @seealso
#' \code{\link{vprofile_windExtinction}}
#' 
#' @examples
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Load example plot plant data
#' data(exampleforestMED)
#' 
#' #Canopy height (in m)
#' h= max(exampleforestMED$treeData$Height/100) 
#' d0 = 0.67*h
#' z0 = 0.08*h
#' 
#' #Height values (cm)
#' z = seq(50,1000, by=50)
#' zm = z/100 # (in m)
#' 
#' # Leaf area density
#' lad = vprofile_leafAreaDensity(exampleforestMED, SpParamsMED, draw = FALSE,
#'                                z = c(0,z))
#'   
#' # Effective drag
#' Cd = 0.2
#' Cx = Cd*lad
#'   
#' # canopy turbulence model
#' wind_canopyTurbulenceModel(zm, Cx,h,d0,z0)
#' 
#' @name wind
wind_canopyTurbulenceModel <- function(zm, Cx, hm, d0, z0, model = "k-epsilon") {
    .Call(`_medfate_windCanopyTurbulenceModel`, zm, Cx, hm, d0, z0, model)
}

#' @rdname wind
wind_canopyTurbulence <- function(zmid, LAD, canopyHeight, u, windMeasurementHeight = 200, model = "k-epsilon") {
    .Call(`_medfate_windCanopyTurbulence`, zmid, LAD, canopyHeight, u, windMeasurementHeight, model)
}

.windSpeedAtCanopyHeight <- function(wind20H, canopyHeight) {
    .Call(`_medfate_windSpeedAtCanopyHeight`, wind20H, canopyHeight)
}

.unshelteredMidflameWindSpeed <- function(wind20H, fuelBedHeight) {
    .Call(`_medfate_unshelteredMidflameWindSpeed`, wind20H, fuelBedHeight)
}

.shelteredMidflameWindSpeed <- function(wind20H, crownFillProportion, topCanopyHeight) {
    .Call(`_medfate_shelteredMidflameWindSpeed`, wind20H, crownFillProportion, topCanopyHeight)
}

#' @rdname fuel_properties
#' 
#' @param topShrubHeight Shrub stratum top height (in m).
#' @param bottomCanopyHeight Canopy base height (in m).
#' @param topCanopyHeight Canopy top height (in m).
#' @param canopyCover Canopy percent cover.
fuel_windAdjustmentFactor <- function(topShrubHeight, bottomCanopyHeight, topCanopyHeight, canopyCover) {
    .Call(`_medfate_windAdjustmentFactor`, topShrubHeight, bottomCanopyHeight, topCanopyHeight, canopyCover)
}

.windSpeedAtHeightOverCanopy <- function(z, wind20H, canopyHeight) {
    .Call(`_medfate_windSpeedAtHeightOverCanopy`, z, wind20H, canopyHeight)
}

.windExtinctionProfile <- function(z, wind20H, LAIc, canopyHeight) {
    .Call(`_medfate_windExtinctionProfile`, z, wind20H, LAIc, canopyHeight)
}

#' Wood formation
#' 
#' Functions to initialize and expand a ring of tracheids to simulate secondary growth.
#' 
#' @param ring An object of class \code{\link{ring}} returned by function \code{woodformation_initRing}.
#' @param psi Water potential (in MPa).
#' @param Tc Temperature in Celsius.
#' @param Nc Number of active cells in the cambium.
#' @param phi0 Initial value of cell extensibility (in MPa-1 day-1)
#' @param pi0 Initial value of cell osmotic potential (in MPa)
#' @param CRD0 Initial value of cell radial diameter
#' @param Y_P Turgor pressure yield threshold (in MPa)
#' @param Y_T Temperature yield threshold (in Celsius)
#' @param h Cell wall hardening coefficient (in day-1)
#' @param s Cell wall softening coefficient (unitless)
#' @param pi Osmotic potential (in MPa)
#' @param phi Cell extensibility (in MPa-1 day-1)
#' @param DHa,DSd,DHd Enthalpy of activation, enthalpy difference and entropy difference between the catalytically active and inactive states of the enzymatic system (Parent et al. 2010).
#'
#' @return 
#' Function \code{woodformation_initRing()} returns a list of class 'ring', 
#' that is a list containing a data frame \code{cells} and two vectors: \code{P} and \code{SA}. 
#' Dataframe \code{cells} contains the columns "formation_date", "phi", "pi" and "CRD" and as many rows as dates processed. 
#' Vectors \code{P} and \code{SA} contain, respectively, the number of cells produced and the sapwood area 
#' corresponding to the ring of cells (assuming a tangencial radius of 20 micrometers). 
#' 
#' Function \code{woodformation_growRing()} modifies the input 'ring' object according to the environmental conditions given as input.
#' 
#' Function \code{woodformation_relativeExpansionRate()} returns a numeric scalar with the relative expansion rate. 
#' 
#' Function \code{woodformation_temperatureEffect()} returns a scalar between 0 and 1 reflecting the temperature effect on tissue formation rate. 
#' 
#' Function \code{woodformation_relativeGrowthRate} returns the annual growth rate, relative to cambium perimeter, estimated from initial and final diameter values.
#' 
#' @note Code modified from package xylomod by Antoine Cabon, available at GitHub
#' 
#' @author
#' Antoine Cabon, CTFC
#' 
#' Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @references
#' Cabon A, \enc{Fernández-de-Uña}{Fernandez-de-Una} L, Gea-Izquierdo G, Meinzer FC, Woodruff DR, \enc{Martínez-Vilalta}{Martinez-Vilalta} J, De \enc{Cáceres}{Caceres} M. 2020a. Water potential control of turgor-driven tracheid enlargement in Scots pine at its xeric distribution edge. New Phytologist 225: 209–221.
#' 
#' Cabon A, Peters RL, Fonti P, \enc{Martínez-Vilalta}{Martinez-Vilalta}  J, De \enc{Cáceres}{Caceres} M. 2020b. Temperature and water potential co-limit stem cambial activity along a steep elevational gradient. New Phytologist: nph.16456.
#' 
#' Parent, B., O. Turc, Y. Gibon, M. Stitt, and F. Tardieu. 2010. Modelling temperature-compensated physiological rates, based on the co-ordination of responses to temperature of developmental processes. Journal of Experimental Botany 61:2057–2069.
#' 
#' @seealso \code{\link{growth}}
#' 
#' @name woodformation
woodformation_initRing <- function() {
    .Call(`_medfate_initialize_ring`)
}

#' @rdname woodformation
woodformation_temperatureEffect <- function(Tc, Y_T = 5.0, DHa = 87.5e3, DSd = 1.09e3, DHd = 333e3) {
    .Call(`_medfate_temperature_function`, Tc, Y_T, DHa, DSd, DHd)
}

#' @rdname woodformation
woodformation_relativeExpansionRate <- function(psi, Tc, pi, phi, Y_P, Y_T) {
    .Call(`_medfate_relative_expansion_rate`, psi, Tc, pi, phi, Y_P, Y_T)
}

#' @rdname woodformation
woodformation_growRing <- function(ring, psi, Tc, Nc = 8.85, phi0 = 0.13, pi0 = -0.8, CRD0 = 8.3, Y_P = 0.05, Y_T = 5.0, h = 0.043*1.8, s = 1.8) {
    invisible(.Call(`_medfate_grow_ring`, ring, psi, Tc, Nc, phi0, pi0, CRD0, Y_P, Y_T, h, s))
}

# Register entry points for exported C++ functions
methods::setLoadAction(function(ns) {
    .Call(`_medfate_RcppExport_registerCCallable`)
})