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

#' @rdname aspwb
#' @keywords internal
aspwbInput <- function(crop_factor, control, soil) {
    .Call(`_medfate_aspwbInput`, crop_factor, control, soil)
}

#' @rdname communication
#' @keywords internal
aspwb_day_inner <- function(internalCommunication, x, date, meteovec, latitude, elevation, slope = NA_real_, aspect = NA_real_, runon = 0.0, lateralFlows = NULL, waterTableDepth = NA_real_, modifyInput = TRUE) {
    .Call(`_medfate_aspwb_day_inner`, internalCommunication, x, date, meteovec, latitude, elevation, slope, aspect, runon, lateralFlows, waterTableDepth, modifyInput)
}

#' @rdname aspwb
#' @keywords internal
aspwb_day <- function(x, date, meteovec, latitude, elevation, slope = NA_real_, aspect = NA_real_, runon = 0.0, lateralFlows = NULL, waterTableDepth = NA_real_, modifyInput = TRUE) {
    .Call(`_medfate_aspwb_day`, x, date, meteovec, latitude, elevation, slope, aspect, runon, lateralFlows, waterTableDepth, modifyInput)
}

.defineASPWBDailyOutput <- function(latitude, elevation, slope, aspect, dateStrings, x) {
    .Call(`_medfate_defineASPWBDailyOutput`, latitude, elevation, slope, aspect, dateStrings, x)
}

.fillASPWBDailyOutput <- function(l, x, sDay, iday) {
    invisible(.Call(`_medfate_fillASPWBDailyOutput`, l, x, sDay, iday))
}

#' Simulation in agricultural areas
#'
#' Function \code{aspwb_day} performs water balance for a single day in an agriculture location.
#' Function \code{aspwb} performs water balance for multiple days in an agriculture location.
#' 
#' @param crop_factor Agriculture crop factor.
#' @param soil An object of class \code{\link{data.frame}} or \code{\link{soil}}.
#' @param control A list with default control parameters (see \code{\link{defaultControl}}).
#' @param x An object of class \code{\link{aspwbInput}}.
#' @param meteo A data frame with daily meteorological data series (see \code{\link{spwb}}). 
#' @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 lateralFlows Lateral source/sink terms for each soil layer (interflow/to from adjacent locations) as mm/day.
#' @param waterTableDepth Water table depth (in mm). When not missing, capillarity rise will be allowed if lower than total soil depth.
#' @param modifyInput Boolean flag to indicate that the input \code{x} object is allowed to be modified during the simulation.
#' @param waterTableDepth Water table depth (in mm). When not missing, capillarity rise will be allowed if lower than total soil depth.
#'   
#' @author
#' Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{spwbInput}}, \code{\link{spwb}}  
#' 
#' @examples
#' 
#' control <- defaultControl()
#' examplesoil <- defaultSoilParams(4)
#' 
#' x <- aspwbInput(0.75, control, examplesoil)
#' 
#' # Day to be simulated
#' d <- 100
#' meteovec <- unlist(examplemeteo[d,-1])
#' date <- as.character(examplemeteo$dates[d])
#' 
#' #Call simulation function for a single days
#' sd <- aspwb_day(x, date, meteovec,  
#'                latitude = 41.82592, elevation = 100) 
#' 
#' #Call simulation function for multiple days
#' S <- aspwb(x, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' @name aspwb
#' @keywords internal
aspwb <- function(x, meteo, latitude, elevation, slope = NA_real_, aspect = NA_real_, waterTableDepth = NA_real_) {
    .Call(`_medfate_aspwb`, x, meteo, latitude, elevation, slope, aspect, waterTableDepth)
}

#' Physical and biophysical utility functions
#' 
#' Internal 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
#' @keywords internal
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).
#' @keywords internal
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).
#' 
#' @keywords internal
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).
#' @keywords internal
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).
#' @keywords internal
biophysics_leafVapourPressure <- function(leafTemp, leafPsi) {
    .Call(`_medfate_leafVapourPressure`, leafTemp, leafPsi)
}

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

#' @rdname biophysics
#' @param temp Temperature (ºC).
#' @keywords internal
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
#' @keywords internal
carbon_sugarStarchDynamicsLeaf <- function(sugarConc, starchConc, eqSugarConc) {
    .Call(`_medfate_sugarStarchDynamicsLeaf`, sugarConc, starchConc, eqSugarConc)
}

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

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

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

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

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

#' @rdname carbon
#' @keywords internal
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
#' @keywords internal
carbon_sapwoodStructuralLivingBiomass <- function(SA, H, L, V, woodDensity, conduit2sapwood) {
    .Call(`_medfate_sapwoodStructuralLivingBiomass`, SA, H, L, V, woodDensity, conduit2sapwood)
}

#' @rdname carbon
#' @keywords internal
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).
#' @keywords internal
carbon_carbonCompartments <- function(x, biomassUnits = "g_m2") {
    .Call(`_medfate_carbonCompartments`, x, biomassUnits)
}

#' Internal communication
#'
#' Functions for internal communication. Not to be called by users.
#' 
#' @param internalCommunication List for internal communication.
#' @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 lateralFlows Lateral source/sink terms for each soil layer (interflow/to from adjacent locations) as mm/day.
#' @param waterTableDepth Water table depth (in mm). When not missing, capillarity rise will be allowed if lower than total soil depth.
#' @param modifyInput Boolean flag to indicate that the input \code{x} object is allowed to be modified during the simulation.
#' @param model String for model, either "spwb" or "growth".
#' 
#' @name communication
#' @keywords internal
copy_model_output <- function(internalCommunication, x, model) {
    .Call(`_medfate_copyModelOutput`, internalCommunication, x, model)
}

#' @rdname communication
#' @keywords internal
general_communication_structures <- function(numCohorts, nlayers, ncanlayers, ntimesteps, model) {
    .Call(`_medfate_generalCommunicationStructures`, numCohorts, nlayers, ncanlayers, ntimesteps, model)
}

#' @rdname communication
#' @keywords internal
instance_communication_structures <- function(x, model) {
    .Call(`_medfate_instanceCommunicationStructures`, x, model)
}

.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(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Calculate fuel properties according to FCCS
#' fccs <- fuel_FCCS(exampleforest, SpParamsMED)
#' fccs
#'   
#' #Calculate fire behavior according to FCCS
#' fire_FCCS(fccs)
#'   
#'  
#' @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
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
fire_leafThermalFactor <- function(SLA, h = 130.0, c = 2500.0) {
    .Call(`_medfate_leafThermalFactor`, SLA, h, c)
}

#' @rdname fire_severity
#' @keywords internal
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
#' @keywords internal
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 fillWithGenus A boolean flag to try imputation of missing values using genus 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(exampleforest)
#' 
#' #A plant-level way to obtain stand basal area
#' sum(plant_basalArea(exampleforest, SpParamsMED), na.rm=TRUE)
#' 
#' #The analogous plant-level function for LAI
#' sum(plant_LAI(exampleforest, SpParamsMED))
#'   
#' #The analogous plant-level function for fuel loading
#' sum(plant_fuelLoading(exampleforest, SpParamsMED))
#'       
#' #Summary function for 'forest' objects can be also used
#' summary(exampleforest, SpParamsMED)
#' 
#' #Cohort IDs in the models
#' plant_ID(exampleforest, SpParamsMED)
#'       
#' @name plant_values
#' @keywords internal
plant_ID <- function(x, SpParams, treeOffset = 0L, shrubOffset = 0L) {
    .Call(`_medfate_cohortIDs`, x, SpParams, treeOffset, shrubOffset)
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

#' @rdname plant_values
#' @param competitionEffect Logical flag to indicate the inclusion of competition effect on LAI estimates.
#' @keywords internal
plant_LAI <- function(x, SpParams, gdd = NA_real_, bounded = TRUE, competitionEffect = TRUE) {
    .Call(`_medfate_cohortLAI`, x, SpParams, gdd, bounded, competitionEffect)
}

#' 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
#' @keywords internal
herb_foliarBiomass <- function(x, SpParams) {
    .Call(`_medfate_herbFoliarBiomass`, x, SpParams)
}

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

#' @rdname herb_values
#' @keywords internal
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 fillWithGenus A boolean flag to try imputation of missing values using genus 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(exampleforest)
#' 
#' # Species basal area in the forest plot
#' species_basalArea(exampleforest, SpParamsMED)
#'   
#' # Value of parameter "Psi_Extract" for two species
#' species_parameter(c("Pinus halepensis", "Quercus ilex"), SpParamsMED, "Psi_Extract")
#'     
#' @name species_values
#' @keywords internal
species_basalArea <- function(x, SpParams) {
    .Call(`_medfate_speciesBasalArea`, x, SpParams)
}

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

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

#' @rdname species_values
#' @keywords internal
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
#' @keywords internal
stand_fuelLoading <- function(x, SpParams, gdd = NA_real_, includeDead = TRUE) {
    .Call(`_medfate_standFuelLoading`, x, SpParams, gdd, includeDead)
}

#' @rdname stand_values
#' @keywords internal
stand_shrubVolume <- function(x, SpParams) {
    .Call(`_medfate_standShrubVolume`, x, SpParams)
}

#' @rdname stand_values
#' @keywords internal
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)
}

#' Input for simulation models (deprecated)
#'
#' Functions \code{forest2spwbInput()} and \code{forest2growthInput()} take an object of class \code{\link{forest}} 
#' and a soil data input to 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{data.frame}} or \code{\link{soil}}, containing soil parameters per soil layer.
#' @param control A list with default control parameters (see \code{\link{defaultControl}}).
#' 
#' @details
#' Function \code{forest2aboveground()} extracts height and species identity from plant cohorts of \code{x}, 
#' and calculate leaf area index and crown ratio. 
#' 
#' \emph{IMPORTANT NOTE}: Function names \code{forest2spwbInput()} and \code{forest2growthInput()} are now internal and deprecated, but 
#' they can still be used for back-compatibility. They correspond to functions \code{\link{spwbInput}} and \code{\link{growthInput}} 
#' 
#' @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}.}
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso \code{\link{spwbInput}}, \code{\link{soil}},  
#' \code{\link{forest}}, \code{\link{SpParamsMED}}, \code{\link{defaultSoilParams}}, \code{\link{plant_ID}}
#' 
#' @examples
#' #Load example plot plant data
#' data(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' # Aboveground parameters
#' forest2aboveground(exampleforest, SpParamsMED)
#' 
#' # Example of aboveground parameters taken from a forest
#' # described using LAI and crown ratio
#' data(exampleforest2)
#' forest2aboveground(exampleforest2, SpParamsMED)
#' 
#' # Define soil with default soil params (4 layers)
#' examplesoil <- defaultSoilParams(4)
#'
#' # Bewowground parameters (distribution of fine roots)
#' forest2belowground(exampleforest, examplesoil, SpParamsMED)
#' 
#' 
#' @name forest2aboveground
#' @keywords internal
forest2aboveground <- function(x, SpParams, gdd = NA_real_, loading = FALSE) {
    .Call(`_medfate_forest2aboveground`, x, SpParams, gdd, loading)
}

#' @rdname forest2aboveground
#' @keywords internal
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). 
#' 
#' @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 loadingOffset A vector of length five with fine fuel loading values (canopy, shrub, herb, woody and litter) to be added to loading estimations from \code{forest}.
#' @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).}
#'   }
#'   
#' 
#' 
#' @references
#' 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(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Show stratification of fuels
#' fuel_stratification(exampleforest, SpParamsMED)
#'   
#' #Calculate fuel properties according to FCCS
#' fccs <- fuel_FCCS(exampleforest, SpParamsMED)
#' fccs
#' 
#' 
#' @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()), loadingOffset = as.numeric( c(0.0, 0.0, 0.0, 0.0, 0.0)), gdd = NA_real_, heightProfileStep = 10.0, maxHeightProfile = 5000, bulkDensityThreshold = 0.05, depthMode = "crownaverage") {
    .Call(`_medfate_FCCSproperties`, object, SpParams, cohortFMC, loadingOffset, 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}}
#' 
#' @keywords internal
mortality_dailyProbability <- function(stressValue, stressThreshold) {
    .Call(`_medfate_dailyMortalityProbability`, stressValue, stressThreshold)
}

#' @rdname communication
#' @keywords internal
growth_day_inner <- function(internalCommunication, x, date, meteovec, latitude, elevation, slope = NA_real_, aspect = NA_real_, runon = 0.0, lateralFlows = NULL, waterTableDepth = NA_real_, modifyInput = TRUE) {
    invisible(.Call(`_medfate_growthDay_inner`, internalCommunication, x, date, meteovec, latitude, elevation, slope, aspect, runon, lateralFlows, waterTableDepth, modifyInput))
}

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

.defineGrowthDailyOutput <- function(latitude, elevation, slope, aspect, dateStrings, x) {
    .Call(`_medfate_defineGrowthDailyOutput`, latitude, elevation, slope, aspect, dateStrings, x)
}

.fillGrowthDailyOutput <- function(l, x, sDay, iday) {
    invisible(.Call(`_medfate_fillGrowthDailyOutput`, l, x, sDay, iday))
}

#' 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}).
#' @param waterTableDepth Water table depth (in mm). When not missing, capillarity rise will be allowed if lower than total soil depth.
#' 
#' @details
#' Detailed model description is available in the medfate book. 
#' Simulations using the 'Sperry' or 'Sureau' 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 = "Sureau"}; 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 = "Sureau"}; 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 = "Sureau"}; 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
#' \donttest{
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#'   
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#'   
#' #Initialize soil with default soil params (4 layers)
#' examplesoil <- defaultSoilParams(4)
#' 
#' #Initialize model input
#' x1 <- growthInput(exampleforest, examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' G1 <- growth(x1, examplemeteo, latitude = 41.82592, elevation = 100)
#'  
#' #Switch to 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' 
#' #Initialize model input
#' x2 <- growthInput(exampleforest,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' G2 <-growth(x2, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' #Switch to 'Sureau' transpiration mode
#' control <- defaultControl("Sureau")
#' 
#' #Initialize model input
#' x3 <- growthInput(exampleforest,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' G3 <-growth(x3, examplemeteo, latitude = 41.82592, elevation = 100)
#' }
#'       
growth <- function(x, meteo, latitude, elevation, slope = NA_real_, aspect = NA_real_, CO2ByYear = numeric(0), waterTableDepth = NA_real_) {
    .Call(`_medfate_growth`, x, meteo, latitude, elevation, slope, aspect, CO2ByYear, waterTableDepth)
}

#' 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(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Initialize soil with default soil params (4 layers)
#' examplesoil <- defaultSoilParams(4)
#' 
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#' 
#' #Switch to 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' 
#' #Initialize input
#' x <- spwbInput(exampleforest,examplesoil, SpParamsMED, control)
#' 
#' #Leaf vulnerability curves
#' hydraulics_vulnerabilityCurvePlot(x, type="leaf")
#' 
#' #Stem vulnerability curves
#' hydraulics_vulnerabilityCurvePlot(x, type="stem")
#'              
#' @name hydraulics_conductancefunctions
#' @keywords internal
hydraulics_psi2K <- function(psi, psi_extract, exp_extract = 3.0) {
    .Call(`_medfate_Psi2K`, psi, psi_extract, exp_extract)
}

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

#' @rdname hydraulics_conductancefunctions
#' @keywords internal
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_xylemConductanceSigmoid <- function(psi, kxylemmax, P50, slope) {
    .Call(`_medfate_xylemConductanceSigmoid`, psi, kxylemmax, P50, slope)
}

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

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

#' @rdname hydraulics_conductancefunctions
#' @keywords internal
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
#' @keywords internal
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 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 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 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 dE Increment of 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 whether negative flux (i.e. from plant to soil) is allowed.
#' @param maxNsteps Maximum number of steps in the construction of supply functions.
#' 
#' @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
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
hydraulics_ECrit <- function(psiUpstream, kxylemmax, c, d, pCrit = 0.001) {
    .Call(`_medfate_ECrit`, psiUpstream, kxylemmax, c, d, pCrit)
}

#' @rdname hydraulics_supplyfunctions
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
hydraulics_E2psiBelowground <- function(E, hydraulicNetwork, psiIni = as.numeric( c(0))) {
    .Call(`_medfate_E2psiBelowground`, E, hydraulicNetwork, psiIni)
}

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

#' @rdname hydraulics_supplyfunctions
#' @keywords internal
hydraulics_E2psiNetwork <- function(E, hydraulicNetwork, psiIni = as.numeric( c(0))) {
    .Call(`_medfate_E2psiNetwork`, E, hydraulicNetwork, psiIni)
}

#' @rdname hydraulics_supplyfunctions
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
hydraulics_supplyFunctionBelowground <- function(hydraulicNetwork, minFlow = 0.0, pCrit = 0.001) {
    .Call(`_medfate_supplyFunctionBelowground`, hydraulicNetwork, minFlow, pCrit)
}

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

#' @rdname hydraulics_supplyfunctions
#' @keywords internal
hydraulics_supplyFunctionNetwork <- function(hydraulicNetwork, minFlow = 0.0, pCrit = 0.001) {
    .Call(`_medfate_supplyFunctionNetwork`, hydraulicNetwork, minFlow, pCrit)
}

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

#' @rdname hydraulics_supplyfunctions
#' @keywords internal
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 root_P50,root_slope Parameters of the Sigmoid function for roots (root xylem vulnerability curve).
#' @param stem_P50,stem_slope Parameters of the Sigmoid function for stems (stem xylem vulnerability curve).
#' @param leaf_P50,leaf_slope Parameters of the Sigmoid 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
#' @keywords internal
hydraulics_maximumSoilPlantConductance <- function(krhizomax, krootmax, kstemmax, kleafmax) {
    .Call(`_medfate_maximumSoilPlantConductance`, krhizomax, krootmax, kstemmax, kleafmax)
}

#' @rdname hydraulics_scalingconductance
#' @keywords internal
hydraulics_soilPlantResistancesSigmoid <- function(psiSoil, psiRhizo, psiStem, PLCstem, psiLeaf, PLCleaf, krhizomax, n, alpha, krootmax, root_P50, root_slope, kstemmax, stem_P50, stem_slope, kleafmax, leaf_P50, leaf_slope) {
    .Call(`_medfate_soilPlantResistancesSigmoid`, psiSoil, psiRhizo, psiStem, PLCstem, psiLeaf, PLCleaf, krhizomax, n, alpha, krootmax, root_P50, root_slope, kstemmax, stem_P50, stem_slope, kleafmax, leaf_P50, leaf_slope)
}

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

#' @rdname hydraulics_scalingconductance
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
hydraulics_taperFactorSavage <- function(height) {
    .Call(`_medfate_taperFactorSavage`, height)
}

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

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

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

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

#' Hydraulic-related defoliation
#' 
#' Functions to calculate the proportion of crown defoliation due to hydraulic disconnection.
#'
#' @param psiLeaf Leaf water potential (in MPa).
#' @param c,d Parameters of the Weibull function.
#' @param P50,slope Parameters of the Sigmoid function.
#' @param PLC_crit Critical leaf PLC corresponding to defoliation
#' @param P50_cv Coefficient of variation (in percent) of leaf P50, to describe the
#' variability in hydraulic vulnerability across crown leaves.
#' 
#' @details The functions assume that crowns are made of a population of leaves whose
#' hydraulic vulnerability (i.e. the water potential corresponding to 50% loss of conductance) 
#' follows a Gaussian distribution centered on the input P50 and with a known coefficient of variation (\code{P50_cv}).
#' The slope parameter (or the c exponent in the case of a Weibull function) is considered constant.
#' Leaves are hydraulically disconnected, and shedded, when their embolism rate exceeds a critical value (\code{PLC_crit}).
#' 
#' @return The proportion of crown defoliation.
#' 
#' @author 
#' Hervé Cochard, INRAE
#' 
#' Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso
#' \code{\link{hydraulics_conductancefunctions}}
#' 
#' @name hydraulics_defoliation
#' @keywords internal
hydraulics_proportionDefoliationSigmoid <- function(psiLeaf, P50, slope, PLC_crit = 0.88, P50_cv = 10.0) {
    .Call(`_medfate_proportionDefoliationSigmoid`, psiLeaf, P50, slope, PLC_crit, P50_cv)
}

#' @name hydraulics_defoliation
#' @keywords internal
hydraulics_proportionDefoliationWeibull <- function(psiLeaf, c, d, PLC_crit = 0.88, P50_cv = 10.0) {
    .Call(`_medfate_proportionDefoliationWeibull`, psiLeaf, c, d, PLC_crit, P50_cv)
}

#' @param month Month of the year (from 1 to 12).
#' @param prec Precipitation for a given day (mm).
#' @param rainfallIntensityPerMonth A vector with twelve positions with average intensity of rainfall (in mm/h) for each month.
#' 
#' @rdname hydrology_interception
#' @keywords internal
hydrology_rainfallIntensity <- function(month, prec, rainfallIntensityPerMonth) {
    .Call(`_medfate_rainfallIntensity`, month, prec, rainfallIntensityPerMonth)
}

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

.hydrology_interceptionLiuDay <- function(Rainfall, Cm, p, ER = 0.05) {
    .Call(`_medfate_interceptionLiuDay`, Rainfall, Cm, p, ER)
}

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

#' Bare soil evaporation and herbaceous transpiration
#'
#' Functions:
#' \itemize{
#'   \item{Function \code{hydrology_soilEvaporationAmount} calculates the amount of evaporation from bare soil, following Ritchie (1972).}
#'   \item{Function \code{hydrology_soilEvaporation} calculates the amount of evaporation from bare soil and distributes it among soil layers.}
#'   \item{Function \code{hydrology_herbaceousTranspiration} calculates the amount of transpiration due to herbaceous plants.}
#' }
#' 
#' @param soil An object of class \code{\link{soil}}.
#' @param snowpack The amount of snow (in water equivalents, mm) in the snow pack.
#' @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.
#' 
#' 
#' @return 
#' Function \code{hydrology_soilEvaporationAmount} returns the amount of water evaporated from the soil. 
#' 
#' Function \code{hydrology_soilEvaporation} returns a vector of water evaporated from each soil layer.
#' 
#' @references 
#' Ritchie (1972). Model for predicting evaporation from a row crop with incomplete cover. - Water resources research.
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso  \code{\link{spwb}}, \code{\link{hydrology_waterInputs}}, \code{\link{hydrology_infiltration}}
#' 
#' 
#' @name hydrology_soilEvaporation
#' @keywords internal
hydrology_soilEvaporation <- function(soil, snowpack, soilFunctions, pet, LgroundSWR, modifySoil = TRUE) {
    .Call(`_medfate_soilEvaporation`, soil, snowpack, soilFunctions, pet, LgroundSWR, modifySoil)
}

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

#' Soil infiltration
#'
#' Soil infiltration functions:
#' \itemize{
#'   \item{Function \code{hydrology_infiltrationBoughton} calculates the amount of water that infiltrates into the topsoil, according to the USDA SCS curve number method (Boughton 1989).}
#'   \item{Function \code{hydrology_infiltrationGreenAmpt} calculates the amount of water that infiltrates into the topsoil, according to the model by Green & Ampt (1911).}
#'   \item{Function \code{hydrology_infiltrationAmount} uses either Green & Ampt (1911) or Boughton (1989) to estimate infiltration.}
#'   \item{Function \code{hydrology_infiltrationRepartition} distributes infiltration among soil layers depending on macroporosity.}
#' }
#' 
#' @param input A numeric vector of (daily) water input (in mm of water).
#' @param Ssoil Soil water storage capacity (can be referred to topsoil) (in mm of water).
#' 
#' 
#' @return 
#' Functions \code{hydrology_infiltrationBoughton}, \code{hydrology_infiltrationGreenAmpt} and \code{hydrology_infiltrationAmount} 
#' return the daily amount of water that infiltrates into the soil (in mm of water). 
#' 
#' Function \code{hydrology_infiltrationRepartition} returns the amount of infiltrated water that reaches each soil layer. 
#' 
#' @references 
#' Boughton (1989). A review of the USDA SCS curve number method. - Australian Journal of Soil Research 27: 511-523.
#' 
#' Green, W.H. and Ampt, G.A. (1911) Studies on Soil Physics, 1: The Flow of Air and Water through Soils. The Journal of Agricultural Science, 4, 1-24. 
#' 
#' @details
#'  When using function \code{hydrology_infiltrationGreenAmpt}, the units of \code{Ksat}, \code{t} and \code{psi_wat} have to be in the same system (e.g. cm/h, h and cm). 
#' 
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso  \code{\link{spwb}}, \code{\link{hydrology_waterInputs}}
#' 
#' @name hydrology_infiltration
#' @keywords internal
hydrology_infiltrationBoughton <- function(input, Ssoil) {
    .Call(`_medfate_infiltrationBoughton`, input, Ssoil)
}

#' @rdname hydrology_infiltration
#' 
#' @param t Time of the infiltration event
#' @param psi_w Matric potential at the wetting front
#' @param Ksat hydraulic conductivity at saturation
#' @param theta_sat volumetric content at saturation
#' @param theta_dry volumetric content at the dry side of the wetting front
#' 
#' @keywords internal
hydrology_infiltrationGreenAmpt <- function(t, psi_w, Ksat, theta_sat, theta_dry) {
    .Call(`_medfate_infitrationGreenAmpt`, t, psi_w, Ksat, theta_sat, theta_dry)
}

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

#' @rdname hydrology_infiltration
#' 
#' @param rainfallInput Water from the rainfall event reaching the soil surface (mm)
#' @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 rainfallIntensity rainfall intensity rate (mm/h)
#' @param model Infiltration model, either "GreenAmpt1911" or "Boughton1989"
#' @param K_correction Correction for saturated conductivity, to account for increased infiltration due to macropore presence
#' 
#' @keywords internal
hydrology_infiltrationAmount <- function(rainfallInput, rainfallIntensity, soil, soilFunctions, model = "GreenAmpt1911", K_correction = 1.0) {
    .Call(`_medfate_infiltrationAmount`, rainfallInput, rainfallIntensity, soil, soilFunctions, model, K_correction)
}

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

#' Water vertical inputs
#' 
#' High-level functions to define water inputs into the soil of a stand:
#'  
#' \itemize{
#'   \item{Function \code{hydrology_waterInputs} performs canopy water interception and snow accumulation/melt.}
#'   \item{Function \code{hydrology_snowMelt} estimates snow melt using a simple energy balance, according to Kergoat (1998).}
#' }
#' 
#' @param x An object of class \code{\link{spwbInput}} or \code{\link{growthInput}}.
#' @param prec Precipitation for the given day (mm)
#' @param pet Potential evapotranspiration for the given day (mm)
#' @param rainfallIntensity Rainfall intensity rate (mm/h).
#' @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-active radiation (PAR) reaching the ground.
#' @param LgroundSWR Percentage of short-wave radiation (SWR) reaching the ground.
#' @param modifyInput Boolean flag to indicate that the input \code{x} object should be modified during the simulation.
#' 
#' @details 
#' The function simulates different vertical hydrological processes, which are described separately in other functions. 
#' If \code{modifyInput = TRUE} the function will modify the \code{x} object (including both soil moisture and 
#' the snowpack on its surface) as a result of simulating hydrological processes.
#' 
#' @return 
#' Function \code{hydrology_waterInputs} 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{Snowmelt}{Snow melted during the day, and added to the water infiltrated.}
#' \item{NetRain}{Rainfall reaching the ground.}
#' 
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @references
#'  Kergoat L. (1998). A model for hydrological equilibrium of leaf area index on a global scale. Journal of Hydrology 212–213: 268–286.
#' 
#' @seealso \code{\link{spwb_day}}, \code{\link{hydrology_rainInterception}}, \code{\link{hydrology_soilEvaporation}}
#' 
#' 
#' @name hydrology_verticalInputs
#' @keywords internal
hydrology_waterInputs <- function(x, prec, rainfallIntensity, pet, tday, rad, elevation, Cm, LgroundPAR, LgroundSWR, modifyInput = TRUE) {
    .Call(`_medfate_waterInputs`, x, prec, rainfallIntensity, pet, tday, rad, elevation, Cm, LgroundPAR, LgroundSWR, modifyInput)
}

#' Soil water balance
#' 
#' Function \code{hydrology_soilWaterBalance} estimates water balance of soil layers given water inputs/outputs, including the simulation of water movement within the soil.
#' 
#' @param soil 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 rainfallInput Amount of water from rainfall event (after excluding interception), in mm.
#' @param rainfallIntensity Rainfall intensity, in mm/h.
#' @param snowmelt Amount of water originated from snow melt, in mm.
#' @param sourceSink Local source/sink term for each soil layer (from soil evaporation or plant transpiration/redistribution)
#'        as mm/day.
#' @param runon Surface water amount running on the target area from upslope (in mm).
#' @param lateralFlows Lateral source/sink terms for each soil layer (interflow/to from adjacent locations) as mm/day.
#' @param waterTableDepth Water table depth (in mm). When not missing, capillarity rise will be allowed if lower than total soil depth.
#' @param infiltrationMode Infiltration model, either "GreenAmpt1911" or "Boughton1989"
#' @param infiltrationCorrection Correction for saturated conductivity, to account for increased infiltration due to macropore presence
#' @param soilDomains Either "buckets" (multi-bucket domain), "single" (for single-domain Richards) or "dual" (for dual-permeability model).
#' @param nsteps Number of time steps per day
#' @param max_nsubsteps Maximum number of substeps per time step
#' @param modifySoil Boolean flag to indicate that the input \code{soil} object should be modified during the simulation.
#' 
#' @seealso  \code{\link{spwb}}, \code{\link{hydrology_waterInputs}}, \code{\link{hydrology_infiltration}}
#' 
#' @details
#' The multi-bucket model adds/substracts water to each layer and if content is above field capacity the excess percolates to the layer below.
#' If there is still an excess for the bottom layer, the model will progressively fill upper layers (generating saturation excess if the first layer 
#' becomes saturated). Every day the some layers are over field capacity, the model simulates deep drainage.
#' 
#' The single-domain model simulates water flows by solving Richards's equation using the predictor-corrector method, as described in 
#' Bonan et al. (2019).
#' 
#' The dual-permeability model is an implementation of the model MACRO 5.0 (Jarvis et al. 1991; Larsbo et al. 2005).
#' 
#' Both the multi-bucket and the single-domain model can apply a correction to the infiltration rate to account for macroporosity in infiltration. In the
#' dual-permeability model extra infiltration through macropores is simulated explicitly.
#' 
#' @author 
#' Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' \enc{María González Sanchís}{Maria Gonzalez Sanchis}, UPV-CTFC
#' 
#' @return
#'   Returns a named vector with different elements, depending on \code{soilDomains}. 
#'   
#'   If \code{soilDomains == "buckets"}:
#'   \itemize{
#'     \item{\code{Snowmelt}: Snowmelt input (mm).}
#'     \item{\code{Source/sinks}: Sum of source/sink input across layers (mm).}
#'     \item{\code{Infiltration}: Water infiltrated into the soil (mm).}
#'     \item{\code{InfiltrationExcess}: Excess infiltration in the topmost layer (mm) leading to an increase in runoff.}
#'     \item{\code{SaturationExcess}: Excess saturation in the topmost layer (mm) leading to an increase in runoff.}
#'     \item{\code{Runoff}: Surface runoff generated by saturation excess or infiltration excess (mm).}
#'     \item{\code{DeepDrainage}: Water draining from the bottom layer (mm). This quantity is corrected to close the water balance.}
#'     \item{\code{CapillarityRise}: Water entering the soil via capillarity rise (mm) from the water table, if \code{waterTableDepth} is supplied.}
#'   }
#'   If \code{soilDomains == "single"} the named vector contains the following additional elements:
#'   \itemize{
#'     \item{\code{Correction}: Amount of water (mm) added to deep drainage to correct the water balance.}
#'     \item{\code{VolumeChange}: Change in soil water volume (mm).}
#'     \item{\code{Substep}: Time step of the moisture solving (seconds).}
#'   }
#'  If \code{soilDomains == "dual"} the named vector contains the following additional elements:
#'   \itemize{
#'     \item{\code{Lateral flows}: Sum of water circulating between micropores and macropores, positive when filling micropores (mm).}
#'     \item{\code{InfiltrationMatrix}: Water infiltrated into the soil matrix (mm).}
#'     \item{\code{InfiltrationMacropores}: Water infiltrated into the soil macropore domain (mm).}
#'     \item{\code{InfiltrationExcessMatrix/InfiltrationExcessMacropores}: Excess infiltration in the topmost layer (mm) leading to an increase in runoff.}
#'     \item{\code{SaturationExcessMatrix/SaturationExcessMacropores}: Excess saturation in the topmost layer (mm) leading to an increase in runoff.}
#'     \item{\code{DrainageMatrix}: Water draining from the bottom layer of the matrix domain (mm). This quantity is corrected to close water balance in the micropore domain.}
#'     \item{\code{DrainageMacropores}: Water draining from the bottom layer of the macropore domain (mm). This quantity is corrected to close the water balance in the macropore domain.}
#'     \item{\code{CorrectionMatrix}: Amount of water (mm) added to deep drainage of soil matrix to correct the water balance.}
#'     \item{\code{CorrectionMacropores}: Amount of water (mm) added to deep drainage of macropores to correct the water balance.}
#'     \item{\code{MatrixVolumeChange}: Change in soil water volume in the soil matrix domain (mm).}
#'     \item{\code{MacroporeVolumeChange}: Change in soil water volume in the macropore domain (mm).}
#'   }
#' 
#' @references
#' Bonan, G. (2019). Climate change and terrestrial ecosystem modeling. Cambridge University Press, Cambridge, UK. 
#' 
#' Jarvis, N.J., Jansson, P‐E., Dik, P.E. & Messing, I. (1991). Modelling water and solute transport in macroporous soil. I. Model description and sensitivity analysis. Journal of Soil Science, 42, 59–70. 
#' 
#' Larsbo, M., Roulier, S., Stenemo, F., Kasteel, R. & Jarvis, N. (2005). An Improved Dual‐Permeability Model of Water Flow and Solute Transport in the Vadose Zone. Vadose Zone Journal, 4, 398–406. 
#' 
#' @examples
#' # Define soil parameters
#' spar <- defaultSoilParams(4)
#' 
#' # Initializes soil hydraulic parameters
#' examplesoil <- soil(spar)
#' 
#' # Water balance in a multi-bucket model
#' hydrology_soilWaterBalance(examplesoil, "VG", 10, 5, 0, c(-1,-1,-1,-1), 
#'                            soilDomains = "buckets", modifySoil = FALSE)
#'                            
#' # Water balance in a single-domain model (Richards equation)
#' hydrology_soilWaterBalance(examplesoil, "VG", 10, 5, 0, c(-1,-1,-1,-1), 
#'                            soilDomains = "single", modifySoil = FALSE)
#'                     
#' # Water balance in a dual-permeability model (MACRO)
#' hydrology_soilWaterBalance(examplesoil, "VG", 10, 5, 0, c(-1,-1,-1,-1), 
#'                            soilDomains = "dual", modifySoil = FALSE)
#'   
#' @name hydrology_soilWaterBalance
#' @keywords internal
hydrology_soilWaterBalance <- function(soil, soilFunctions, rainfallInput, rainfallIntensity, snowmelt, sourceSink, runon = 0.0, lateralFlows = NULL, waterTableDepth = NA_real_, infiltrationMode = "GreenAmpt1911", infiltrationCorrection = 5.0, soilDomains = "buckets", nsteps = 24L, max_nsubsteps = 3600L, modifySoil = TRUE) {
    .Call(`_medfate_soilWaterBalance`, soil, soilFunctions, rainfallInput, rainfallIntensity, snowmelt, sourceSink, runon, lateralFlows, waterTableDepth, infiltrationMode, infiltrationCorrection, soilDomains, nsteps, max_nsubsteps, modifySoil)
}

.gammln <- function(xx) {
    .Call(`_medfate_gammln`, xx)
}

.betacf <- function(a, b, x) {
    .Call(`_medfate_betacf`, a, b, x)
}

.incbeta <- function(a, b, x) {
    .Call(`_medfate_incbeta`, a, b, x)
}

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

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

#' @rdname hydraulics_supplyfunctions
#' @keywords internal
hydraulics_initSperryNetworks <- function(x) {
    .Call(`_medfate_initSperryNetworks`, x)
}

#' 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(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Define soil with default soil params (4 layers)
#' examplesoil <- defaultSoilParams(4)
#' 
#' #Initialize control parameters
#' control <- defaultControl(transpirationMode="Sperry")
#' 
#' #Initialize input
#' x2 <- spwbInput(exampleforest,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)
}

#' Sureau-ECOS inner functions for testing only
#' 
#' Function \code{initSureauNetworks} 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 = "Sureau"}.
#'  
#' @return Function \code{initSureauNetworks} 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
#' @keywords internal
initSureauNetworks <- function(x) {
    .Call(`_medfate_initSureauNetworks`, x)
}

#' @rdname sureau_ecos
#' @param network A hydraulic network element of the list returned by \code{initSureauNetworks}
#' @param dt Smallest time step (seconds)
#' @param opt Option flag vector
#' @param stemCavitationRecovery,leafCavitationRecovery A string indicating how refilling of embolized conduits is done:
#'           \itemize{
#'             \item{"none" - no refilling.}
#'             \item{"annual" - every first day of the year.}
#'             \item{"rate" - following a rate of new sapwood formation.}
#'             \item{"total" - instantaneous complete refilling.}
#'           }
#' @keywords internal
semi_implicit_integration <- function(network, dt, opt, stemCavitationRecovery = "annual", leafCavitationRecovery = "total") {
    invisible(.Call(`_medfate_semi_implicit_integration`, network, dt, opt, stemCavitationRecovery, leafCavitationRecovery))
}

#' Advanced radiation transfer functions
#' 
#' Functions \code{light_layerDirectIrradianceFraction} and \code{light_layerDiffuseIrradianceFraction} calculate 
#' the fraction of above-canopy direct and diffuse radiation 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 absorbed 
#' 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 ZF Fraction of sky angles.
#' @param Ib0 Above-canopy direct incident radiation.
#' @param Id0 Above-canopy diffuse incident radiation.
#' @param leafAngle Average leaf inclination angle (in radians).
#' @param leafAngleSD Standard deviation of leaf inclination angle (in radians).
#' @param p,q Parameters of the beta distribution for leaf angles
#' @param ClumpingIndex The extent to which foliage has a nonrandom spatial distribution.
#' @param alpha A vector of leaf absorbance by species.
#' @param gamma A vector of leaf reflectance values.
#' @param solarElevation Solar elevation (in radians).
#' @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).
#' 
#' @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_layerDirectIrradianceFraction}, \code{light_layerDiffuseIrradianceFraction}
#' 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}}, \code{\link{light_basic}}
#' 
#' @examples
#' solarElevation <- 0.67 # in radians
#' SWR_direct <- 1100
#' SWR_diffuse <- 300
#' PAR_direct <- 550
#' PAR_diffuse <- 150
#' 
#' LAI <- 2
#' nlayer <- 10
#' LAIlayerlive <- matrix(rep(LAI/nlayer,nlayer),nlayer,1)
#' LAIlayerdead <- matrix(0,nlayer,1)
#' meanLeafAngle <- 60 # in degrees
#' sdLeafAngle <- 20
#' 
#' beta <- light_leafAngleBetaParameters(meanLeafAngle*(pi/180), sdLeafAngle*(pi/180))
#' 
#' ## Extinction coefficients
#' kb <- light_directionalExtinctionCoefficient(beta["p"], beta["q"], solarElevation)
#' kd_PAR <- 0.5
#' kd_SWR <- kd_PAR/1.35
#' @name light_advanced
#' @keywords internal
light_leafAngleCDF <- function(leafAngle, p, q) {
    .Call(`_medfate_leafAngleCDF`, leafAngle, p, q)
}

#' @rdname light_advanced
#' @keywords internal
light_leafAngleBetaParameters <- function(leafAngle, leafAngleSD) {
    .Call(`_medfate_leafAngleBetaParameters`, leafAngle, leafAngleSD)
}

#' @rdname light_advanced
#' @keywords internal
light_directionalExtinctionCoefficient <- function(p, q, solarElevation) {
    .Call(`_medfate_directionalExtinctionCoefficient`, p, q, solarElevation)
}

#' @rdname light_advanced
#' @keywords internal
light_layerDirectIrradianceFraction <- function(LAIme, LAImd, LAImx, kb, ClumpingIndex, alpha, gamma, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_layerDirectIrradianceFraction`, LAIme, LAImd, LAImx, kb, ClumpingIndex, alpha, gamma, trunkExtinctionFraction)
}

#' @rdname light_advanced
#' @keywords internal
light_layerDiffuseIrradianceFraction <- function(LAIme, LAImd, LAImx, K, ClumpingIndex, ZF, alpha, gamma, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_layerDiffuseIrradianceFraction`, LAIme, LAImd, LAImx, K, ClumpingIndex, ZF, alpha, gamma, trunkExtinctionFraction)
}

#' @rdname light_advanced
#' @keywords internal
light_cohortSunlitShadeAbsorbedRadiation <- function(Ib0, Id0, LAIme, LAImd, LAImx, kb, K, ClumpingIndex, ZF, alpha, gamma, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_cohortSunlitShadeAbsorbedRadiation`, Ib0, Id0, LAIme, LAImd, LAImx, kb, K, ClumpingIndex, ZF, alpha, gamma, trunkExtinctionFraction)
}

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

#' @rdname light_advanced
#' @keywords internal
light_instantaneousLightExtinctionAbsortion <- function(LAIme, LAImd, LAImx, p, q, ClumpingIndex, alphaSWR, gammaSWR, ddd, ntimesteps = 24L, trunkExtinctionFraction = 0.1) {
    .Call(`_medfate_instantaneousLightExtinctionAbsortion`, LAIme, LAImd, LAImx, p, q, ClumpingIndex, alphaSWR, gammaSWR, ddd, ntimesteps, trunkExtinctionFraction)
}

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

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

#' Radiation extinction functions used in basic transpiration sub-model
#' 
#' @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.
#' 
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @seealso  \code{\link{spwb}}, \code{\link{light_advanced}}
#' 
#' @name light_basic
#' @keywords internal
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_basic
#' @keywords internal
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_basic
#' @keywords internal
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_basic
#' @keywords internal
light_cohortAbsorbedSWRFraction <- function(z, x, SpParams, gdd = NA_real_) {
    .Call(`_medfate_cohortAbsorbedSWRFraction`, z, x, SpParams, gdd)
}

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

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

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

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

.rootDistributionComplete <- function(x, SpParams, fillMissingRootParams) {
    .Call(`_medfate_rootDistributionComplete`, x, SpParams, fillMissingRootParams)
}

#' Input for simulation models
#'
#' Functions \code{spwbInput()} and \code{growthInput()} take an object of class \code{\link{forest}} 
#' and a soil data input to create input objects for simulation functions \code{\link{spwb}} (or \code{\link{pwb}}) and \code{\link{growth}}, respectively. 
#' 
#' @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 soil An object of class \code{\link{data.frame}} or \code{\link{soil}}, containing soil parameters per soil layer.
#' @param control A list with default control parameters (see \code{\link{defaultControl}}).
#' 
#' @details
#' Functions \code{spwbInput()} and \code{growthInput()} initialize inputs differently depending on control parameters.
#' 
#' \emph{IMPORTANT NOTE}: Older function names \code{\link{forest2spwbInput}} and \code{\link{forest2growthInput}} are now deprecated, but 
#' they can still be used for back-compatibility.
#' 
#' @return 
#' Function \code{spwbInput()} 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{soil}: A data frame with initialized soil parameters (see \code{\link{soil}}).}
#'     \item{\code{snowpack}: The amount of snow (in mm) in the snow pack over the soil.}
#'     \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 = "Sureau"} 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 = "Sureau"} 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 (Weibull).}
#'         \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_co2}: Coefficient regulating the influence of atmospheric CO2 concentration on gross photosynthesis.}
#'         \item{\code{WUE_vpd}: 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-1·MPa-1).}
#'         \item{\code{Kmax_rootxylem}: Sapwood-specific hydraulic conductivity of root xylem (in kg H2O·s-1·m-1·MPa-1).}
#'         \item{\code{VCleaf_kmax}: Maximum leaf hydraulic conductance (in mmol H2O·s-1·m-2·MPa-1).}
#'         \item{\code{VCleaf_c}, \code{VCleaf_d}: Parameters of the leaf vulnerability curve (Weibull).}
#'         \item{\code{VCstem_kmax}: Maximum stem xylem conductance (in mmol H2O·s-1·m-2·MPa-1).}
#'         \item{\code{VCstem_c}, \code{VCstem_d}: Parameters of the stem xylem vulnerability curve (Weibull).}
#'         \item{\code{VCroot_c}, \code{VCroot_d}: Parameters of the root xylem vulnerability curve (Weibull).}
#'         \item{\code{Plant_kmax}: Maximum whole-plant conductance (in mmol H2O·s-1·m-2·MPa-1).}
#'         \item{\code{FR_leaf}, \code{FR_stem}, \code{FR_root}: Fraction of whole-plant resistance corresponding to each segment.}
#'       }
#'       If \code{control$transpirationMode = "Sureau"} 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{Gsw_AC_slope}: Slope of the Gsw vs Ac/Cs relationship (see \code{\link{photo_photosynthesisBaldocchi}}).}
#'         \item{\code{Gs_P50}, \code{Gs_slope}: Parameters of the curve describing the decrease in stomatal conductance as a function of leaf water potential (sigmoid).}
#'         \item{\code{Vmax298}: Maximum Rubisco carboxylation 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-1·MPa-1).}
#'         \item{\code{Kmax_rootxylem}: Sapwood-specific hydraulic conductivity of root xylem (in kg H2O·s-1·m-1·MPa-1).}
#'         \item{\code{VCleaf_kmax}: Maximum leaf hydraulic conductance (in mmol H2O·s-1·m-2·MPa-1).}
#'         \item{\code{VCleaf_c}, \code{VCleaf_d}: Parameters of the leaf vulnerability curve (Weibull).}
#'         \item{\code{VCleaf_P50}, \code{VCleaf_slope}: Parameters of the leaf vulnerability curve (sigmoid).}
#'         \item{\code{VCstem_kmax}: Maximum stem xylem conductance (in mmol H2O·s-1·m-2·MPa-1).}
#'         \item{\code{VCstem_c}, \code{VCstem_d}: Parameters of the stem xylem vulnerability curve (Weibull).}
#'         \item{\code{VCstem_P50}, \code{VCstem_slope}: Parameters of the stem xylem vulnerability curve (sigmoid).}
#'         \item{\code{VCroot_c}, \code{VCroot_d}: Parameters of the root xylem vulnerability curve (Weibull).}
#'         \item{\code{VCroot_P50}, \code{VCroot_slope}: Parameters of the root xylem vulnerability curve (sigmoid).}
#'         \item{\code{Plant_kmax}: Maximum whole-plant conductance (in mmol H2O·s-1·m-2·MPa-1).}
#'         \item{\code{FR_leaf}, \code{FR_stem}, \code{FR_root}: Fraction of whole-plant resistance corresponding to each segment.}
#'       }
#'     }
#'     \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{growthInput()} 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(exampleforest)
#' 
#' # Example of aboveground parameters taken from a forest
#' # described using LAI and crown ratio
#' data(exampleforest2)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' 
#' # Define soil with default soil params (4 layers)
#' examplesoil <- defaultSoilParams(4)
#' 
#' # Initialize control parameters using 'Granier' transpiration mode
#' control <- defaultControl("Granier")
#' 
#' # Prepare spwb input
#' spwbInput(exampleforest, examplesoil, SpParamsMED, control)
#'                 
#' # Prepare input for 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' spwbInput(exampleforest,examplesoil,SpParamsMED, control)
#' 
#' # Prepare input for 'Sureau' transpiration mode
#' control <- defaultControl("Sureau")
#' spwbInput(exampleforest,examplesoil,SpParamsMED, control)
#' 
#' # Example of initialization from a forest 
#' # described using LAI and crown ratio
#' control <- defaultControl("Granier")
#' spwbInput(exampleforest2, examplesoil, SpParamsMED, control)
#' 
#' @name modelInput
spwbInput <- function(x, soil, SpParams, control) {
    .Call(`_medfate_spwbInput`, x, soil, SpParams, control)
}

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

#' @rdname forest2aboveground
#' @keywords internal
forest2spwbInput <- function(x, soil, SpParams, control) {
    .Call(`_medfate_forest2spwbInput`, x, soil, SpParams, control)
}

#' @rdname forest2aboveground
#' @keywords internal
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, fillWithGenus = TRUE) {
    .Call(`_medfate_speciesNumericParameterWithImputation`, species, SpParams, parName, fillMissing, fillWithGenus)
}

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

.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
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
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
#' 
#' @keywords internal
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
#' @keywords internal
photo_electronLimitedPhotosynthesis <- function(Q, Ci, GT, Jmax) {
    .Call(`_medfate_electronLimitedPhotosynthesis`, Q, Ci, GT, Jmax)
}

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

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

#' @rdname photo
#' @param Gsw_AC_slope Slope of the An/C vs Gsw relationship 
#' @param Gsw_AC_intercept Intercept of the An/C vs Gsw relationship 
#' @keywords internal
photo_photosynthesisBaldocchi <- function(Q, Catm, Tleaf, u, Vmax298, Jmax298, leafWidth, Gsw_AC_slope, Gsw_AC_intercept) {
    .Call(`_medfate_photosynthesisBaldocchi`, Q, Catm, Tleaf, u, Vmax298, Jmax298, leafWidth, Gsw_AC_slope, Gsw_AC_intercept)
}

#' @rdname photo
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
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 Z100 A vector of depths (in mm) corresponding to 100% 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{spwbInput}}, \code{\link{soil}}
#'
#' @examples
#' #Load example plot plant data
#' data(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' ntree <- nrow(exampleforest$treeData)
#' 
#' #Initialize soil with default soil params
#' s <- defaultSoilParams(4)
#' 
#' #Calculate conic root system for trees
#' V1 <- root_conicDistribution(Z=rep(2000,ntree), s$widths)            
#' print(V1)
#'      
#' #Calculate LDR root system for trees (Schenck & Jackson 2002)
#' V2 <- root_ldrDistribution(Z50 = rep(200,ntree), 
#'                            Z95 = rep(1000,ntree),
#'                            Z100 = rep(NA, ntree), s$widths)
#' print(V2)     
#' 
#' @name root
#' @keywords internal
root_conicDistribution <- function(Zcone, d) {
    .Call(`_medfate_conicDistribution`, Zcone, d)
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

#' @rdname root
#' @keywords internal
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 Initialized 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_saturatedWaterDepth()} returns the depth to saturation 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
#' 
#' # Define soil with default params
#' soil_df <- defaultSoilParams(4)
#' soil_df
#' 
#' # Initialize soil parameters and state variables
#' s = soil(soil_df)
#' 
#' # Plot Saxton's and Van Genuchten's water retention curves
#' soil_retentionCurvePlot(s, model="both")
#' 
#' @name soil_texture
#' @keywords internal
soil_saturatedConductivitySX <- function(clay, sand, bd, om = NA_real_, mmol = TRUE) {
    .Call(`_medfate_saturatedConductivitySaxton`, clay, sand, bd, om, mmol)
}

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

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

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

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

#' @rdname soil_texture
#' @param ksat saturated hydraulic conductance
#' @keywords internal
soil_psi2kVG <- function(ksat, n, alpha, theta_res, theta_sat, psi) {
    .Call(`_medfate_psi2kVanGenuchten`, ksat, n, alpha, theta_res, theta_sat, psi)
}

#' @rdname soil_texture
#' @keywords internal
soil_psi2cVG <- function(n, alpha, theta_res, theta_sat, psi) {
    .Call(`_medfate_psi2cVanGenuchten`, n, alpha, theta_res, theta_sat, psi)
}

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

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

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

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

#' @rdname soil_texture
#' @keywords internal
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
#' @keywords internal
soil_waterSAT <- function(soil, model = "SX") {
    .Call(`_medfate_waterSAT`, soil, model)
}

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

#' @rdname soil_texture
#' @keywords internal
soil_waterPsi <- function(soil, psi, model = "SX") {
    .Call(`_medfate_waterPsi`, soil, psi, model)
}

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

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

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

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

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

#' @rdname soil_texture
#' @keywords internal
soil_conductivity <- function(soil, model = "SX") {
    .Call(`_medfate_conductivity`, soil, model)
}

#' @rdname soil_texture
#' @keywords internal
soil_capacitance <- function(soil, model = "SX") {
    .Call(`_medfate_capacitance`, soil, model)
}

#' @rdname soil_texture
#' @keywords internal
soil_saturatedWaterDepth <- function(soil, model = "SX") {
    .Call(`_medfate_saturatedWaterDepth`, soil, model)
}

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

#' @rdname soil_texture
#' @keywords internal
soil_campbellParamsClappHornberger <- function(soilType) {
    .Call(`_medfate_campbellParamsClappHornberger`, soilType)
}

#' @rdname soil_texture
#' @keywords internal
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 x 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).
#' 
#' @return
#' Function \code{soil} returns a data frame of class \code{soil} with the following columns:
#' \itemize{
#'   \item{\code{widths}: 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{nitrogen}: Sum of total nitrogen (ammonia, organic and reduced nitrogen) for each layer (in g/kg).}
#'   \item{\code{rfc}: Percentage of rock fragment content for each layer.}
#'   \item{\code{macro}: Macroporosity for each layer (estimated using Stolf et al. 2011).}
#'   \item{\code{Ksat}: Saturated soil conductivity for each layer (estimated using function \code{\link{soil_saturatedConductivitySX}}.}
#'   \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{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.}
#' }
#' 
#' @author Miquel De \enc{Cáceres}{Caceres} Ainsa, CREAF
#' 
#' @details 
#' Function \code{summary} 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. 
#' 
#' If available, the user can specify columns \code{VG_alpha}, \code{VG_n}, \code{VG_theta_res}, \code{VG_theta_sat} and \code{K_sat},
#' to override Van Genuchten parameters an saturated conductivity estimated from pedotransfer functions when calling function \code{soil}. 
#' 
#' @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_redefineLayers}}, \code{\link{soil_psi2thetaSX}}, \code{\link{soil_psi2thetaVG}}, \code{\link{spwb}}, \code{\link{defaultSoilParams}}
#' 
#' @examples
#' # Default parameters
#' df_soil <- defaultSoilParams()
#' 
#' # Initializes soil
#' s = soil(df_soil)
#' s
#' 
#' # Prints soil characteristics according to Saxton's water retention curve
#' summary(s, model="SX")
#' 
#' # Prints soil characteristics according to Van Genuchten's water retention curve
#' summary(s, model="VG")
#' 
#' # Add columns 'VG_theta_sat' and 'VG_theta_res' with custom values
#' df_soil$VG_theta_sat <- 0.400 
#' df_soil$VG_theta_res <- 0.040 
#' 
#' # Reinitialize soil (should override estimations)
#' s2 = soil(df_soil)
#' s2
#' summary(s2, model="VG")
#' @name soil
soil <- function(x, VG_PTF = "Toth") {
    .Call(`_medfate_soilInit`, x, VG_PTF)
}

.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 widths 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_SAT Relative water content (in percent volume) at saturation 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).
#' @param tstep Time step (interval) in seconds.
#' 
#' @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
#' #Define soil and complete parameters
#' examplesoil = soil(defaultSoilParams(4))
#' 
#' 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
#' @keywords internal
soil_thermalCapacity <- function(soil, model = "SX") {
    .Call(`_medfate_thermalCapacity`, soil, model)
}

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

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

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

.getWeatherDates <- function(meteo) {
    .Call(`_medfate_getWeatherDates`, meteo)
}

.defineSPWBDailyOutput <- function(latitude, elevation, slope, aspect, dateStrings, x) {
    .Call(`_medfate_defineSPWBDailyOutput`, latitude, elevation, slope, aspect, dateStrings, x)
}

.fillSPWBDailyOutput <- function(l, x, sDay, iday) {
    invisible(.Call(`_medfate_fillSPWBDailyOutput`, l, x, sDay, iday))
}

#' 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 called \code{"dates"} or \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}).
#' @param waterTableDepth Water table depth (in mm). When not missing, capillarity rise will be allowed if lower than total soil depth.
#' 
#' @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 'Sureau' transpiration mode is illustrated by function \code{\link{transp_transpirationSureau}} and was described for model SurEau-Ecos v2.0 in Ruffault et al. (2022).} 
#' }
#' Simulations using the 'Sperry' or 'Sureau' 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 rainfall (in mm).}
#'     \item{\code{"Snow"}: Precipitation as snowfall (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{"InfiltrationExcess"}: Excess infiltration in the topmost layer leading to an increase in runoff (in mm).}
#'     \item{\code{"SaturationExcess"}: Excess saturation in the topmost layer leading to an increase in runoff (in mm).}
#'     \item{\code{"CapillarityRise"}: Water entering the soil via capillarity rise (mm) from the water table, if \code{waterTableDepth} is supplied.}
#'     \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 = "Sureau"}).}
#'   \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 = "Sureau"}).}
#'   \item{\code{"Soil"}: A list with the following subelements:}
#'   \itemize{
#'     \item{\code{"SWC"}: Soil water content (percent of soil volume) in each soil layer (and overall).}
#'     \item{\code{"RWC"}: Relative soil moisture content (relative to field capacity) in each soil layer (and overall).}
#'     \item{\code{"REW"}: Relative extractable water (min. psi = -5 MPa) in each soil layer (and overall).}
#'     \item{\code{"ML"}: Soil water volume in each soil layer (in L/m2) (and overall).}
#'     \item{\code{"Psi"}: Soil water potential in each soil layer (in MPa) (and overall).}
#'     \item{\code{"PlantExt"}: Plant extraction from each soil layer (in mm) (and overall).}
#'     \item{\code{"HydraulicInput"}: Water that entered the layer coming from other layers and transported via the plant hydraulic network (in mm) (and overall).}
#'   }
#'   \item{\code{"Snow"}: A data frame where the following variable (in columns) is given for each simulated day (in rows):}
#'   \itemize{
#'     \item{\code{"SWE"}: Snow water equivalent (mm) of the snow pack.}
#'   }
#'   \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{"LeafPLC"}: A data frame with the average daily proportion of leaf conductance loss of each plant (\[0-1\]).}
#'     \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="Sureau"}, 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{"LeafPLC"}: A data frame with the average daily proportion of leaf conductance loss of each plant (\[0-1\]).}
#'     \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{extract}}, \code{\link{summary.spwb}},  \code{\link{forest}}, \code{\link{aspwb}}
#' 
#' @examples
#' \donttest{
#' #Load example daily meteorological data
#' data(examplemeteo)
#' 
#' #Load example plot plant data
#' data(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Define soil with default soil params (4 layers)
#' examplesoil <- defaultSoilParams(4)
#' 
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#' 
#' #Initialize input
#' x1 <- spwbInput(exampleforest,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' S1 <- spwb(x1, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' #Switch to 'Sperry' transpiration mode
#' control <- defaultControl("Sperry")
#' 
#' #Initialize input
#' x2 <- spwbInput(exampleforest,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' S2 <- spwb(x2, examplemeteo, latitude = 41.82592, elevation = 100)
#' 
#' #Switch to 'Sureau' transpiration mode
#' control <- defaultControl("Sureau")
#' 
#' #Initialize input
#' x3 <- spwbInput(exampleforest,examplesoil, SpParamsMED, control)
#' 
#' #Call simulation function
#' S3 <- spwb(x3, examplemeteo, latitude = 41.82592, elevation = 100)
#' }
#'                 
#' @name spwb
spwb <- function(x, meteo, latitude, elevation, slope = NA_real_, aspect = NA_real_, CO2ByYear = numeric(0), waterTableDepth = NA_real_) {
    .Call(`_medfate_spwb`, x, meteo, latitude, elevation, slope, aspect, CO2ByYear, waterTableDepth)
}

#' @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, 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)
}

#' @rdname communication
#' @keywords internal
spwb_day_inner <- function(internalCommunication, x, date, meteovec, latitude, elevation, slope = NA_real_, aspect = NA_real_, runon = 0.0, lateralFlows = NULL, waterTableDepth = NA_real_, modifyInput = TRUE) {
    invisible(.Call(`_medfate_spwbDay_inner`, internalCommunication, x, date, meteovec, latitude, elevation, slope, aspect, runon, lateralFlows, waterTableDepth, modifyInput))
}

#' 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 lateralFlows Lateral source/sink terms for each soil layer (interflow/to from adjacent locations) as mm/day.
#' @param waterTableDepth Water table depth (in mm). When not missing, capillarity rise will be allowed if lower than total soil depth.
#' @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 'Sureau' transpiration mode is illustrated by function \code{\link{transp_transpirationSureau}} and was described for model SurEau-Ecos v2.0 in Ruffault et al. (2022).} 
#' }
#' 
#' Simulations using the 'Sperry' or 'Sureau' 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{"Psi"}: Soil water potential (in MPa) at the end of the day.}
#'       \item{\code{"HerbTranspiration"}: Water extracted by herbaceous plants from each soil layer (in mm).}
#'       \item{\code{"HydraulicInput"}: Water entering each soil layer from other layers, transported via plant roots (in mm).}
#'       \item{\code{"HydraulicOutput"}: Water leaving each soil layer (going to other layers or the transpiration stream) (in mm).}
#'       \item{\code{"PlantExtraction"}: Water extracted by woody plants from each soil layer (in mm).}
#'     }
#'   }
#'   \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 = "Sureau"}:
#' \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(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Define soil parameters
#' examplesoil <- defaultSoilParams(4)
#' 
#' # Day to be simulated
#' d <- 100
#' meteovec <- unlist(examplemeteo[d,-1])
#' date <- as.character(examplemeteo$dates[d])
#' 
#' #Simulate water balance one day only (Granier mode)
#' control <- defaultControl("Granier")
#' x1 <- spwbInput(exampleforest,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 <- spwbInput(exampleforest, examplesoil, SpParamsMED, control)
#' sd2 <-spwb_day(x2, date, meteovec,
#'               latitude = 41.82592, elevation = 100, slope=0, aspect=0)
#' 
#' #Simulate water balance for one day only (Sureau mode)
#' control <- defaultControl("Sureau")
#' x3 <- spwbInput(exampleforest, 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  <- growthInput(exampleforest,examplesoil, SpParamsMED, control)
#' sd4 <- growth_day(x4, date, meteovec,
#'                 latitude = 41.82592, elevation = 100, slope=0, aspect=0)
#' 
#' #Simulate water and carbon balance for one day only (Sperry mode)
#' control <- defaultControl("Sperry")
#' x5  <- growthInput(exampleforest,examplesoil, SpParamsMED, control)
#' sd5 <- growth_day(x5, date, meteovec,
#'                 latitude = 41.82592, elevation = 100, slope=0, aspect=0)
#' 
#' #Simulate water and carbon balance for one day only (Sureau mode)
#' control <- defaultControl("Sureau")
#' x6  <- growthInput(exampleforest,examplesoil, SpParamsMED, control)
#' sd6 <- growth_day(x6, date, meteovec,
#'                 latitude = 41.82592, elevation = 100, slope=0, aspect=0)
#' 
#' @name spwb_day
spwb_day <- function(x, date, meteovec, latitude, elevation, slope = NA_real_, aspect = NA_real_, runon = 0.0, lateralFlows = NULL, waterTableDepth = NA_real_, modifyInput = TRUE) {
    .Call(`_medfate_spwbDay`, x, date, meteovec, latitude, elevation, slope, aspect, runon, lateralFlows, waterTableDepth, modifyInput)
}

#' 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
#' @keywords internal
moisture_sapwoodWaterCapacity <- function(Al2As, height, V, L, wd) {
    .Call(`_medfate_sapwoodWaterCapacity`, Al2As, height, V, L, wd)
}

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

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

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

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

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

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

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

#' @rdname moisture
#' @keywords internal
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.
#' 
#' @keywords internal
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
#' @keywords internal
transp_transpirationSureau <- function(x, meteo, day, latitude, elevation, slope, aspect, canopyEvaporation = 0.0, snowMelt = 0.0, soilEvaporation = 0.0, herbTranspiration = 0.0, modifyInput = TRUE) {
    .Call(`_medfate_transpirationSureau`, 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_transpirationSureau} 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 'Sureau' 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_transpirationSureau}, 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_transpirationSureau}:
#'   \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(exampleforest)
#' 
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Define soil with default soil params (4 layers)
#' examplesoil <- defaultSoilParams(4)
#' 
#' #Initialize control parameters
#' control <- defaultControl("Granier")
#' 
#' #Initialize input
#' x1 <- spwbInput(exampleforest,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 <- spwbInput(exampleforest,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 'Sureau' transpiration mode
#' control <- defaultControl("Sureau")
#' 
#' #Initialize input
#' x3 <- spwbInput(exampleforest,examplesoil, SpParamsMED, control)
#' 
#' # Transpiration according to Sureau model
#' t3 <- transp_transpirationSureau(x3, examplemeteo, 1, 
#'                                   latitude = 41.82592, elevation = 100, slope = 0, aspect = 0,
#'                                   modifyInput = FALSE)
#'                                 
#' @name transp_modes
#' @keywords internal
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(exampleforest)
#' 
#' #Canopy height (in m)
#' h= max(exampleforest$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(exampleforest, 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
#' @keywords internal
wind_canopyTurbulenceModel <- function(zm, Cx, hm, d0, z0, model = "k-epsilon") {
    .Call(`_medfate_windCanopyTurbulenceModel`, zm, Cx, hm, d0, z0, model)
}

#' @rdname wind
#' @keywords internal
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)
}

#' Wind adjustment factor for Rothermel's model
#' 
#' Function fuel_windAdjustmentFactor determines the adjustment factor of wind for surface fires, according to Andrews (2012).
#' 
#' @name fuel_windAdjustmentFactor
#' 
#' @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.
#' 
#' @returns A scalar 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.
#' 
#' @keywords internal
#' 
#' @examples
#' #Load example plot plant data
#'  data(exampleforest)
#'   
#' #Default species parameterization
#' data(SpParamsMED)
#' 
#' #Calculate fuel properties according to FCCS
#' fccs <- fuel_FCCS(exampleforest, SpParamsMED)
#' 
#' # Estimate wind adjustment factor
#' fuel_windAdjustmentFactor(fccs$htc[2], fccs$hbc[1], fccs$htc[1], fccs$cover[1])
#' 
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
#' @keywords internal
woodformation_initRing <- function() {
    .Call(`_medfate_initialize_ring`)
}

#' @rdname woodformation
#' @keywords internal
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
#' @keywords internal
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
#' @keywords internal
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`)
})