# # TODO: maybe move this into a module? on soil discretization?
"""
discretize_soil(;
    Δz_m::Vector{T},
    Rootden_::Function,
    uAux_PSIM_init_kPa::Function,
    u_delta18O_init_permil::Function = ((Δz_m) -> missing),
    u_delta2H_init_permil::Function = ((Δz_m) -> missing))

Manually generate a soil- and root-discretization for LWFBrook90.jl. This function can be
used as alternative to loading an input file `soil_discretization.csv`.

Reqired arguments are a vector with thickness of the discretization cells/layers in meter
`Δ_m`, and functions `Rootden_(Δ_m)` that generates the root density, a function
`uAux_PSIM_init_kPa(Δ_m)` that generates the initial values of ψ (kPa)
for all cells as a function of `Δ_m`.

Optional arguments are functions initial conditions of the two isotopic
signatures based on `Δ_m`.
"""
function discretize_soil(;
    Δz_m::Vector{T},
    Rootden_::Function,
    uAux_PSIM_init_kPa::Function,
    u_delta18O_init_permil::Function = ((Δz_m) -> missing),
    u_delta2H_init_permil::Function = ((Δz_m) -> missing)) where {T<:Number}

    # Derive cell interfaces based on the cell spacing Δz_m
    z_cell_interfaces = -[0; cumsum(Δz_m)]
    z_cell_interfaces = round.(z_cell_interfaces; digits=6)

    # Output required DataFrame
    soil_discretizationDF = DataFrame(
        Upper_m = z_cell_interfaces[Not(end)],
        Lower_m = z_cell_interfaces[Not(1)],
        Rootden_ = Rootden_(Δz_m),
        uAux_PSIM_init_kPa = uAux_PSIM_init_kPa(Δz_m),
        u_delta18O_init_permil = u_delta18O_init_permil(Δz_m),
        u_delta2H_init_permil = u_delta2H_init_permil(Δz_m))

    # TODO: refine Δz if needed for certain parameters
    #       Don't do this if it is a manual definition
    #       But do it if it stems from discretize_soil(::SPAC)

    # Output DataFrame
    return soil_discretizationDF
end

# function setup(soil_horizons::DataFrame, )
#     # checks on structure
#     @assert names(soil_horizons) == ["HorizonNr", "Upper_m", "Lower_m", "shp"]
#     coltypes = [eltype(col) for col in eachcol(input_soil_horizons[1])]
#     @assert coltypes[1] <: Integer
#     @assert coltypes[2] <: Real
#     @assert coltypes[3] <: Real
#     @assert coltypes[4] <: LWFBrook90.AbstractSoilHydraulicParams

#     # return DataFrame(Upper_m, Lower_m, Rootden_, ICs, shp)
#     # where ICs is an Array{State}, and shp an Array{AbstractSoilHydraulicParams}
#     # TODO(bernhard): define
#     return DiscretizedSoilDomain(Upper_m, Lower_m, Rootden_, ICs, shp)
# end

"""
    extend_lowest_horizon(soil_horizons, soil_discretizationDF)

Extends the physical description of the soil domain further downward, it this is
necessary by the requested discretization.
"""
function extend_lowest_horizon(soil_horizons, soil_discretizationDF)
    # Copy to modify
    extended_soil_horizons = deepcopy(soil_horizons)

    # Check extent
    soil_horizon_extent = (maximum(extended_soil_horizons.Upper_m), minimum(extended_soil_horizons.Lower_m))

    @assert soil_horizon_extent[1] ≈ soil_discretizationDF.Upper_m[1] "Spatial domain of soil discretization is not compatible with provided soil horizons"
    if soil_horizon_extent[2] ≈ soil_discretizationDF.Lower_m[end]
        # all okay no need to extend z horizon
    elseif soil_horizon_extent[2] < soil_discretizationDF.Lower_m[end]
        @warn    "Spatial domain soil discretization is smaller than the provided soil horizon information ($soil_horizon_extent m)" * "\n" *
        "         (lower end of requested soil discr.: $(soil_discretizationDF.Lower_m[end]))"
    elseif soil_horizon_extent[2] > soil_discretizationDF.Lower_m[end]
        @warn """
        Spatial domain soil discretization is larger than the provided soil horizon information. Lowest soil horizon will be extended.

        Lowest soil horizon from $(extended_soil_horizons[end,:Upper_m])m to $(extended_soil_horizons[end,:Lower_m])m will be extended down to $(soil_discretizationDF.Lower_m[end])m.
        """
        push!(extended_soil_horizons, extended_soil_horizons[end,:])
        extended_soil_horizons[end,:Lower_m] = soil_discretizationDF.Lower_m[end]
        extended_soil_horizons[end,:Upper_m] = extended_soil_horizons[end-1,:Lower_m]
        extended_soil_horizons[end,:HorizonNr] = extended_soil_horizons[end-1,:HorizonNr] + 1
    end

    return extended_soil_horizons
end

"""
    refine_soil_discretization(
        Δz,
        soil_output_depths_m,
        IDEPTH_m, #=modified_SPAC.pars.params[:IDEPTH_m],
        QDEPTH_m)

Densify discretization of soil domain whenever an interface or an additional cell is needed.
This is the case for interfaces
- when a new soil horizons begins or
- when a state variable needs to be extracted at a specified output depth (add an additional computational cell (consisting of upper AND lower interfaces))
- when a parameter such as IDEPTH_m or QDEPTH_m requires an interface (cell boundary)
"""
function refine_soil_discretization(
    # Δz::Vector,
    prior_soil_discretization::DataFrame, #Δz, TODO: simplify this function so that it uses only Δz, not the full DF soil_discretization
    input_soil_horizons::DataFrame,
    soil_output_depths_m::Vector,
    IDEPTH_m::Real, # = modified_SPAC.pars.params[:IDEPTH_m],
    QDEPTH_m::Real; # = modified_SPAC.pars.params[:QDEPTH_m])

    # Define what is close enough for soil_output_depths_m (also minimal thickness of layers to create with this procedure)
    # ε = 0.005     # thickness of layer to be inserted, [m]
    # ε = 0.010   # thickness of layer to be inserted, [m]
    # ε = 0.020   # thickness of layer to be inserted, [m]
    # ε = 0.025   # thickness of layer to be inserted, [m]
    ε = 0.050   # thickness of layer to be inserted, [m]
    )

    # input_soil_horizons: Physical description of soil domain, with horizons
    # soil_output_depths_m:  Requested depths where the user wants to have the state variables, e.g. to compare with measurements
    # IDEPTH_m depth over which infiltration is distributed, [m] "IDEPTH_m determines the
    #        number of soil layers over which infiltration is distributed when INFEXP is
    #        greater than 0."
    # QDEPTH_m soil depth for SRFL calculation, 0 to prevent SRFL, [m] "QDEPTH_m determines the
    #        layers over which wetness is calculated to determine source area (SRFL)
    #        parameters."



    # IF WE WANT TO SUPPORT Δz as argument instead of soil_discretization_DF:
    # if prior_soil_discretization isa DataFrame
    #     all_lower_interfaces = prior_soil_discretization.Lower_m
    #     # TO DEVELOP: check if this is approximately equal with
    #     # all_lower_interfaces ≈ cumsum(Δz)
    # else
    #     cumsum(Δz)
    #     error("Soil discretization is not a DataFrame. Please, check with the developer why you see this message.")
    # end

    ############
    # 1) Check if discretization needs to be refined
    # 1a) Find out which layers need to be added
    @assert all(abs.(diff(soil_output_depths_m)) .>= ε) """
        Requested soil_output_depths_m (additional layers) must be further away from each other than $ε m.
        Requested were: $(soil_output_depths_m)
    """

    # To accomodate requested soil_output_depths_m
    layers_to_insert = soil_output_depths_m
    interfaces_for_layers = reverse(sort(vcat(layers_to_insert, layers_to_insert .- ε)))

    # To accomodate current parameter settings (e.g. infiltration depth)
    interfaces_for_parameters = sort([
         # interfaces for change in horizons or inflow depths or QDEPTH_m
         input_soil_horizons[!,"Upper_m"];
         input_soil_horizons[!,"Lower_m"];
         -IDEPTH_m; # m (positive) to m (negative)
         -QDEPTH_m;  # m (positive) to m (negative)
         ];
        rev = true)

    all_potentially_needed_interfaces =
        unique(sort([interfaces_for_layers; interfaces_for_parameters]; rev=true))

    # 1b) Find out which interfaces will be too close (< ε)
    existing_interfaces = [0; prior_soil_discretization[!,"Lower_m"]]

    needed_interfaces = zeros(Float64, 0)
    for layer in all_potentially_needed_interfaces
        if !any(abs.(layer .- existing_interfaces) .< ε * 0.5) # 0.5 as
            # insert the layer if none (!any() of the already existing layers is close enough)
            append!(needed_interfaces, layer)
        end
    end
    # Check again that needed_interfaces themselves are spaced widely enough:
    if length(needed_interfaces)>0
        needed_interfaces = needed_interfaces[
            [true;                                    # keep first one for sure
            abs.(diff(needed_interfaces)) .>= ε*0.9]  # keep 2:N layers only if diff large enough (use 0.9ε instead of ε because of floating point precisison)
        ]
    end

   # Add them to the DataFrame and copy the properties of the cell above.
    refined_soil_discr = deepcopy(prior_soil_discretization) # otherwise input argument is modified in-place
    added_interfaces = Float64[] # to report
    for interface_to_add in needed_interfaces
        for i in 1:nrow(refined_soil_discr)
            condition = interface_to_add > refined_soil_discr[i, "Lower_m"]
            #println("$i, interface_to_add:$interface_to_add, $(refined_soil_discr[i, "Lower_m"]) $condition")
            if (condition)
                interface_to_add
                pos = i+1
                insert!.(eachcol(refined_soil_discr),
                        pos,
                        Matrix(refined_soil_discr)[pos-1,:]) # use properties from cell above
                refined_soil_discr[pos-1, "Lower_m"] = interface_to_add[1] # modify
                refined_soil_discr[pos,   "Upper_m"] = interface_to_add[1] # modify
                push!(added_interfaces, interface_to_add[1])
                break
            elseif (i == nrow(refined_soil_discr))
                # If interface to add is below the lowest layer: add it only in case it is about ε lower
                if (interface_to_add >= refined_soil_discr[i, "Lower_m"] - ε)
                    pos = i+1
                    insert!.(eachcol(refined_soil_discr),
                        pos,
                        Matrix(refined_soil_discr)[pos-1,:]) # use properties from cell above
                    refined_soil_discr[pos, "Upper_m"] = refined_soil_discr[pos-1, "Lower_m"]
                    refined_soil_discr[pos, "Lower_m"] = interface_to_add[1]
                    push!(added_interfaces, interface_to_add[1])
                    break
                end
            end
        end
    end
    if length(added_interfaces) != 0
        @warn "Adding interfaces to soil discretization at depths $added_interfaces (m), to allow for either requested IDEPTH/QDEPTH:  $(IDEPTH_m)m/$(QDEPTH_m)m, "*
              "or to accomodate physically defined soil horizons: $(unique(vcat(input_soil_horizons[!,"Upper_m"], input_soil_horizons[!,"Lower_m"]))) (m)"*
              ifelse(length(soil_output_depths_m)==0,".",", or `soil_output_depths_m`=$soil_output_depths_m.")
    end

    # Define IDEPTH_idx and QDEPTH_idx (to be used internally instead of IDEPTH_m and QDEPTH_m)
    is_infiltration_layer_BOOLEAN = -IDEPTH_m .<= refined_soil_discr[!,"Lower_m"]
    is_SRFL_layer_BOOLEAN         = -QDEPTH_m .<= refined_soil_discr[!,"Lower_m"]
    IDEPTH_idx = sum(is_infiltration_layer_BOOLEAN) # lowest node where Lower_m is below IDEPTH_m
    QDEPTH_idx = sum(is_SRFL_layer_BOOLEAN)         # lowest node where Lower_m is below QDEPTH_m

    if (-IDEPTH_m < refined_soil_discr[end,"Lower_m"]) ||
        (-QDEPTH_m < refined_soil_discr[end,"Lower_m"])
        error("QDEPTH_m or IDEPTH_m were defined deeper than the lowest simulation element.")
    end
    ############

    return refined_soil_discr, IDEPTH_idx, QDEPTH_idx
    # return refined_Δz, IDEPTH_idx, QDEPTH_idx
end

"""
    map_soil_horizons_to_discretization(soil_horizons, soil_discretization)

Attribute soil parameters based on the defined horizons to computational cells.
"""
function map_soil_horizons_to_discretization(soil_horizons, computational_grid)

    # Assert expectations
    @assert soil_horizons[1,"Upper_m"] ≈ 0
    @assert computational_grid[1,"Upper_m"] ≈ 0
    @assert soil_horizons[1,"Lower_m"] < 0
    @assert computational_grid[1,"Lower_m"] < 0

    # Find for each soil_discretization node in which horizon it lies
   which_horizon = fill(0, nrow(computational_grid))
    for i = 1:nrow(computational_grid)
        # For each soil discretizations (i) check to which which horizon it belongs
        idx = findfirst(computational_grid[i,"Lower_m"] .>= soil_horizons[:,"Lower_m"])

        if (isnothing(idx) && i == nrow(computational_grid))
            idx = nrow(soil_horizons)
        end

        which_horizon[i] = soil_horizons[idx,"HorizonNr"]
    end
    computational_grid2 = deepcopy(computational_grid) # prevents modification of argument
    computational_grid2[:,"HorizonNr"] = which_horizon

    # # If lowest discretized layer is only ε lower than defined soil horizons extrapolate
    # # the properties of the lowest soil horizon to that infinitesimal thin layer
    # # NOTE: this is not needed anymore as we have the function `extend_lowest_horizon()`
    # if ( (computational_grid2[end,"HorizonNr"] == 0) &
    #     (abs(computational_grid2[end,"Lower_m"] - soil_horizons[end,"Lower_m"]) < ε) )
    #     computational_grid2[end,"HorizonNr"] = nrow(soil_horizons)
    # end

    nlayer_before_join = nrow(computational_grid2)
    computational_grid2 = innerjoin(computational_grid2[:, [:Upper_m, :Lower_m, :Rootden_, :uAux_PSIM_init_kPa, :u_delta18O_init_permil, :u_delta2H_init_permil, :HorizonNr]], # remove old horizon limits and shp
                                    rename(soil_horizons, :Upper_m => :Horizon_Upper_m, :Lower_m => :Horizon_Lower_m, ), # append new horizons and shp
                                    makeunique=false,
                                    #select(soil_horizons, Not([:Upper_m,:Lower_m])),
                                    on = :HorizonNr)
                                    # NOTE: not using leftjoin becaus it transforms types to Union{Missing, Float64} instead of only Float64
                                    # Therefore, use innerjoin and check manually that no rows in computational_grid2 are lost
    @assert nrow(computational_grid2) == nlayer_before_join """
        When merging soil properties defined in input CSV-file containing soil horizons, some layers from
        the soil discretization were lost. This is most likely due to the soil horizons not covering the entire discretization domain.
        (Note that the discretized domain is defined in a separate input CSV-file.)

        Please check and correct the input files.
        """
        # (Note that the discretized domain is defined in a separate input CSV-file and is further extended to include internal variables IDEPTH_m and QDEPTH_m.)
    ############

    # Some final checks and return
    THICK_m        = computational_grid2[!,"Upper_m"] - computational_grid2[!,"Lower_m"] # thickness of soil layer [m]
    THICK          = 1000*(THICK_m)                                 # thickness of soil layer [mm]
    @assert all(THICK_m .> 0.0) "All discretized cells msut be bigger than zero...."
    PSIM_init      = computational_grid2[!,"uAux_PSIM_init_kPa"]    # initial condition PSIM [kPa]
    @assert all(PSIM_init .<= 0) "Initial matrix psi must be negative or zero"

    return computational_grid2
    # TODO: replace in define_LWFB90_p:
    #       - soil_discr["SHP"] with soil_horizons_Dict["soil_discretization"].shp
    #         i.e. parametrizedSPAC....something.shp
    #       - soil_horizons_Dict["NLAYER"] with nrow(soil_horizons_Dict["soil_discretization"])
    #       - soil_horizons_Dict["THICK"] with
    #         i.e. parametrizedSPAC...something.Δz_m*1000
    #       - soil_horizons_Dict["PSIM_init"]
    #       - soil_horizons_Dict["final_Rootden_"]
    #       - soil_horizons_Dict["PSIM_init"]      = soil_horizons_Dict["soil_discretization"][!,"uAux_PSIM_init_kPa"]     # initial condition PSIM [kPa]
    #       - soil_horizons_Dict["d18O_soil_init"] = soil_horizons_Dict["soil_discretization"][!,"u_delta18O_init_permil"] # initial condition soil water δ18O [‰]
    #       - soil_horizons_Dict["d2H_soil_init"]  = soil_horizons_Dict["soil_discretization"][!,"u_delta2H_init_permil"]  # initial condition soil water δ2H [‰]
    #       - soil_horizons_Dict["final_Rootden_"] = soil_horizons_Dict["soil_discretization"][!,"Rootden_"]


    ############
    HEAT = 0 # flag for heat balance; not implemented; parametrizedSPAC.params[:HEAT], @ hardcoded
    # # BELOW IS ONLY UNUSED CODE:... linked to discretization of thermal properties
    # for i = 1:NLAYER
    #     if (HEAT != 0)
    #         error("HEAT must be set to zero, as heat balance is not implemented.")
    #     end
    #     # if (HEAT == 1)
    #     #     # TemperatureNew(i)       = soil_discretization[i,7] we don't have it in the input file!!!
    #     #     if i > NLAYER
    #     #         MUE[i] = THICK[i] / ( THICK[i] + THICK(I+1) )
    #     #         ZL[i]  = 0.5 * ( THICK[i] + THICK(I+1) )
    #     #     else
    #     #         MUE[i] = 0.5
    #     #         ZL[i]  = THICK[i]
    #     #     end
    #     #     TMean[i]   = 0.
    #     # end
    # end

    # # heat flow -------
    # # nmat   = nrow(input_soil_horizons)
    # if (HEAT != 0) error("HEAT must be set to zero, as heat balance is not implemented.") end
    # #       if (HEAT == 1)
    # #        READ (12,*) Comment
    # #        READ (12,*) tTop, Comment
    # #        tTop=tTop
    # #        READ (12,*) tBot, Comment
    # #        tBot=tBot
    # #        READ (12,*) TopInfT, Comment
    # #        READ (12,*) BotInfT, Comment
    # #        READ (12,*) kTopT, Comment
    # #        READ (12,*) kBotT, Comment
    # #        DO 207 I = 1, 7
    # #         READ (12,*) Comment
    # # 207    CONTINUE
    # #        DO 208 I = 1, nmat
    # #         READ (12,*) ilay, SV[i], OV[i], HB1[i], HB2[i], HB3[i]
    # #         TPar[1,I) = SV[i]
    # #         TPar[2,I) = OV[i]
    # #         TPar[3,I) = THDis
    # # C        thermal conductivities -- transfer from [J m-1 s-1 K-1] to  [J mm-1 d-1 K-1]
    # #         TPar[4,I) = HB1[i] * 86.400
    # #         TPar[5,I) = HB2[i] * 86.400
    # #         TPar[6,I) = HB3[i] * 86.400
    # # C         volumetric heat capacities for solid, water and organic -- transfer from [MJ m-2 mm-1 K-1] to [J mm-3 K-1]
    # #         TPar[7,I) = p_CVSOL   # To define in module_CONSTANTS.jl: p_CVSOL = ?  # CVSOL  - volumetric heat capacity of solid (MJ m-2 mm-1 K-1)
    # #         TPar[8,I) = p_CVORG   # To define in module_CONSTANTS.jl: p_CVORG = ?  # volumetric heat capacity of organic material (MJ m-2 mm-1 K-1) (hillel98)
    # #         TPar[9,I) = LWFBrook90.CONSTANTS.p_CVLQ
    # # 208    CONTINUE
    # #        READ (12,*) C
    # #       end

    # # mat       = input_soil_horizons[!,"mat"]
    # # ### # from LWFBrook90R:md_brook90.f95 (lines 105)
    # # TPar = fill(NaN, (nmat, 10))
    # # TPar .= -1
    # # HeatCapOld = fill(NaN, NLAYER)
    # # for i = 1:NLAYER
    # #     HeatCapOld[i] = TPar[mat[i],7] * TPar[mat[i],1] + TPar[mat[i],8]
    # #                 # TPar[2,mat[i])+TPar[9,mat[i])*SWATI[i]/THICK[i]
    # # end
    # ###
end

function overwrite_IC!(soil_discretization_DF, _to_use_IC_soil, simulate_isotopes)
    @assert (keys(_to_use_IC_soil) == (:PSIM_init_kPa, :delta18O_init_permil, :delta2H_init_permil)) | (keys(_to_use_IC_soil) == (:PSIM_init_kPa))
    @assert !((simulate_isotopes) && (keys(_to_use_IC_soil) == (:PSIM_init_kPa)))
    if (keys(_to_use_IC_soil) == (:PSIM_init_kPa, :delta18O_init_permil, :delta2H_init_permil))
        soil_discretization_DF.uAux_PSIM_init_kPa     .= _to_use_IC_soil.PSIM_init_kPa
        soil_discretization_DF.u_delta18O_init_permil .= _to_use_IC_soil.delta18O_init_permil
        soil_discretization_DF.u_delta2H_init_permil  .= _to_use_IC_soil.delta2H_init_permil
    else
        soil_discretization_DF.uAux_PSIM_init_kPa     .= _to_use_IC_soil.PSIM_init_kPa
    end
    return nothing
end

function overwrite_rootden!(soil_discretization_DF, _to_use_root_distribution, _to_use_Δz_thickness_m)
    # remove empty arguments (i.e. nothing)
    _to_use_root_distribution = NamedTuple([k => v for (k,v) in pairs(_to_use_root_distribution) if !isnothing(v)])

    # only possible combinations of arguments
    @assert ((sort(collect(keys(_to_use_root_distribution))) == [:beta, :z_rootMax_m]) |
             (sort(collect(keys(_to_use_root_distribution))) == [:beta]) |
             (sort(collect(keys(_to_use_root_distribution))) == [:root_k, :root_θ_cm, :z_rootMax_m]) |
             (sort(collect(keys(_to_use_root_distribution))) == [:root_k, :root_θ_cm])
    ) "Need to provide complete set of root parameters, i.e. either beta, (beta, z_rootMax_m), (root_k, root_θ_cm, z_rootMax_m), or (root_k, root_θ_cm)"

    # _to_use_root_distribution = (; root_k = 1.0, root_θ_cm = 20); _to_use_Δz_thickness_m = [0.02, 0.02, 0.06, fill(0.1,10)...]
    # _to_use_root_distribution = (; root_k = 1.0); _to_use_Δz_thickness_m = [0.02, 0.02, 0.06, fill(0.1,10)...]
    # _to_use_root_distribution = (; beta = 0.97); _to_use_Δz_thickness_m = [0.02, 0.02, 0.06, fill(0.1,10)...]

    soil_discretization_DF.Rootden_ = LWFBrook90.Rootden_(; _to_use_root_distribution..., Δz_m = _to_use_Δz_thickness_m)
    return nothing
end

"""
    Rootden_(
        [[beta],[root_k, root_θ_cm]];
        Δz_m,
        z_rootMax_m = maximum(cumsum(Δz_m)),
        z_Upper_m = 0)

Define the relative root density in each discretized soil layer either as either
- exponential distribution with mean b = -1cm/ln(β), with input argument `beta` (Gale and Grigal, 1987),
- gamma distribution with shape `root_k`>0 and scale `root_θ`, resulting in a distribution with mean k*θ

This function returns the instantaneous root fraction (dY/dd) (units of -/cm). The values are
    normalized so that the area under the curve sums up to 1.0 (when plotted vs cm).
    The method reduces the amount of roots in a discretization layer in case the effective maximal rooting
    depth comes to lie within that layer - it does not need to modify the discretization in that case.

    (Normalization extends the possible range of β to values larger than 1, and thereby allows
    to describe profiles that are denser at the bottom. Although it is recommended to use the gamma distribution for these cases.)
"""
function Rootden_(;
    beta  = nothing,                       # beta distribution Gale&Grigal 1987
    root_k = nothing, root_θ_cm = nothing, # exponential or gamma distributions
    Δz_m,
    z_rootMax_m = -sum(Δz_m), # total_effective_rooting_depth_m
    z_Upper_m = 0)

    # Notation:
    # z: absolute vertical position (positive upward, i.e. soil values are negative)
    # d: depth below beginning of mineral soil
    @assert z_Upper_m <= 0 "Discretization of humus layer (>0m) is untested."
    @assert minimum(Δz_m) > 0 "Discretization cell height Δz_m must be positive"
    @assert z_rootMax_m < 0 "z_rootMax_m must be negative"

    if !isnothing(beta)
        # Case 1: using β for root parametrization
        @assert isnothing(root_k) && isnothing(root_θ_cm) "beta provided: $beta. You cannot also provided root_k $root_k or root_θ_cm $root_θ_cm."
        # (Gale and Grigal, 1987)
        Y_fct(β, d_cm) = 1 .- β .^ d_cm
        # dYdd_fct(β, d_cm) = -β .^ d_cm .* log(β) # unused
        Y_fct_eval = (d_cm) -> 1 .- beta .^ d_cm

        @assert abs(beta - 1.0) > 0.000001 "Root distribution parameter β is too close to 1.0."*
            "Note that it can be below (or above) 1.0, but not exactly 1. Set it at least 0.000001 away from 1.0"

        # # (Gale and Grigal, 1987) figures 1 and 2
        # d_midpoints_m = -z_Upper_m + Δz_m[1]/2 .+
        #     cumsum(Δz_m) .- Δz_m[1] # d is soil depth d in centimeters
        # # TODO: by using d_midpoint of we make a small error.
        # #       Better would be to use integral between cell boundaries.
        # d_cm = d_midpoints_m*100
        # pl_1 = plot(
        #     Y_fct(0.92, d_cm), d_cm,
        #     xflip = true,
        #     yflip = true, ylabel = "Depth_cm", xlabel = "Cumulative Root Fraction (Y)");
        # plot!(Y_fct(0.95, d_cm), d_cm);
        # plot!(Y_fct(0.97, d_cm), d_cm);
        # pl_2 = plot(
        #     dYdd_fct(0.92, d_cm), d_cm,
        #     yflip = true, ylabel = "Depth_cm", xlabel = "Instantaneous Root Fraction (dY/dd)");
        # plot!(dYdd_fct(0.95, d_cm), d_cm);
        # plot!(dYdd_fct(0.97, d_cm), d_cm);
        # plot(pl_1, pl_2; layout=(1,2))
    elseif !isnothing(root_θ_cm)
        # Case 2: Gamma distribution: using θ (and k) for root parametrization
        @assert isnothing(beta) "root_θ_cm provided: $root_θ_cm. You cannot also provided beta $beta."

        # Gamma distribution: https://en.wikipedia.org/wiki/Gamma_distribution
        # Y_fct(root_k, root_θ_cm, d_cm, z_rootMax_m) = cdf.(truncated.(Gamma.(root_k, root_θ_cm), 0, -100*z_rootMax_m), d_cm)
        Y_fct_eval = (d_cm) -> cdf.(truncated.(Gamma.(root_k, root_θ_cm), 0, -100*z_rootMax_m), d_cm)
    else
        error("No known parameters provided as arguments to `Rootden_()`. Please provide either arguments `beta` or then both: `root_k` and `root_θ_cm`")
    end


    z_interfaces_m = z_Upper_m .+ [0; cumsum(-Δz_m)]
    # 0) Preparatory steps
    # Depth of cell interfaces at which we evaluate the cumulative root fraction
    d_interfaces_cm = -z_interfaces_m * 100
    # # Crop cumulative density at total effective rooting depth
    d_root_depth_cm = -z_rootMax_m * 100

    # 1) Compute cumulative density inside of a cell (i.e. increase in cumulative denisity from upper to lower interface)
    rootden_increase_per_cell = diff(Y_fct_eval(d_interfaces_cm))

    # 2) Correct for total effective rooting depth
    #   a) total rooting depth is below lowest discretization cell                    (nothing needs changing)
    #   b) total rooting depth coincides with lower boundary of a discretization cell (remove roots further below)
    #   c) total rooting depth is inside of a discretization cell                     (remove roots below and correct cell)
    d_upper_interfaces_cm = d_interfaces_cm[Not(end),]
    d_lower_interfaces_cm = d_interfaces_cm[Not(1),]

    # remove roots from cells that are completely below rooting depth
    rootden_increase_per_cell[d_upper_interfaces_cm.>=d_root_depth_cm] .= 0
    rootden_increase_per_cell

    # correct root density in cell where total rooting depth lies
    idx = findfirst(d_root_depth_cm .<= d_lower_interfaces_cm)
    if (!isnothing(idx))
        # First computation above was:
        # Y_fct(β, d_lower_interfaces_cm[idx]) - Y_fct(β, d_upper_interfaces_cm[idx])
        # Corrected computation is:
        #rootden_increase_per_cell[idx] = Y_fct(β, d_root_depth_cm) - Y_fct(β, d_upper_interfaces_cm[idx])
        rootden_increase_per_cell[idx] = Y_fct_eval(d_root_depth_cm) - Y_fct_eval(d_upper_interfaces_cm[idx])
        # (note if d_lower_interfaces_cm[idx] == d_root_depth_cm, this keeps the value before)
    end

    # 3) Normalize the increase in root density over the cell by the cell height
    #    note that this uses Δz_m also for the one cell which contains the total rooting depth
    rootden_increase_per_cmDepth = rootden_increase_per_cell ./ (Δz_m .* 100)

    # 4) Normalize all values so that it sums up to 1.0
    normalizeAUC = (x) -> sum(x) == 0 ? x : x ./ sum(x)
    return normalizeAUC(rootden_increase_per_cmDepth)
end
function Rootden_beta_(
    β;
    Δz_m,
    z_rootMax_m = -sum(Δz_m), # total_effective_rooting_depth_m
    z_Upper_m = 0) # # Upper interface of topmost discretization cell
    # e.g. Δz_m = fill(0.10, 10)
    #      z_rootMax_m = -0.22
    #      β = 0.97
    Rootden_(beta = β, Δz_m = Δz_m, z_rootMax_m = z_rootMax_m, z_Upper_m = z_Upper_m)
end
# Δz_m = fill(0.10, 10) # (m)
# Δz_m = diff(collect(0:0.01:1)) # (m)
# z_rootMax_m = -0.225
# f = Figure(); ax = Axis(f[1, 1], ylabel = "Depth (cm)", xlabel = "Instantaneous root fraction (dY/dd) (Area of 1.0)", yreversed = true)
# Makie.lines!(ax,
#     LWFBrook90.Rootden_(root_β = 0.97, Δz_m = Δz_m), cumsum(Δz_m * 100), label = "No max Root")
# # Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 0.97, Δz_m = Δz_m, z_rootMax_m = z_rootMax_m),
# #     cumsum(Δz_m * 100), label = "maxRoot: $(z_rootMax_m)m")
# # Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 0.97, Δz_m = Δz_m, z_rootMax_m = -1.2),
# #     cumsum(Δz_m * 100), label = "maxRoot: $(-1.2)m")
# # Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 0.97, Δz_m = Δz_m, z_rootMax_m = -0.5),
# #     cumsum(Δz_m * 100), label = "maxRoot: $(-0.5)m",)
# Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 0.97, Δz_m = Δz_m, z_rootMax_m = z_rootMax_m),
#     cumsum(Δz_m * 100), label = "maxRoot: $(z_rootMax_m)m")
# Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 0.97, Δz_m = Δz_m, z_rootMax_m = -0.9),
#     cumsum(Δz_m * 100), label = "maxRoot: $(-0.9)m")
# Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 0.97, Δz_m = Δz_m, z_rootMax_m = -0.5),
#     cumsum(Δz_m * 100), label = "maxRoot: $(-0.5)m",)
# Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 0.999, Δz_m = Δz_m, z_rootMax_m = -0.9),
#     cumsum(Δz_m * 100), label = "β = 0.999")
# # Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 1.0, Δz_m = Δz_m, z_rootMax_m = -0.9),
# #     cumsum(Δz_m * 100), label = "β = 1.0")
# # Makie.lines!(ax, LWFBrook90.Rootden_(root_β = 1.001, Δz_m = Δz_m, z_rootMax_m = -0.9),
# #     cumsum(Δz_m * 100), label = "β = 1.001")
# # Makie.lines!(ax, LWFBrook90.Rootden_(root_k = 1.0, root_θ_cm = -1/log(0.97),
# #     Δz_m = Δz_m, z_rootMax_m = -0.9), cumsum(Δz_m * 100), label = "β = 0.97, k=1.0")
# Makie.lines!(ax, LWFBrook90.Rootden_(root_k = 3.0, root_θ_cm = -1/log(0.97),
#     Δz_m = Δz_m, z_rootMax_m = -0.9), cumsum(Δz_m * 100), label = "β = 0.97, k=3.0")
# Makie.lines!(ax, LWFBrook90.Rootden_(root_k = 5.0, root_θ_cm = -1/log(0.97),
#     Δz_m = Δz_m, z_rootMax_m = -0.9), cumsum(Δz_m * 100), label = "β = 0.97, k=5.0")
# axislegend(ax)
# f



"""
    HammelKennel_transient_root_density()

Take root density distribution parameters and generates root density in time (t) and space (z, negative downward).
"""
function HammelKennel_transient_root_density(;
        timepoints,
        AGE_at_timepoints,
        p_INITRDEP,   p_INITRLEN,
        p_RGROPER_y,  p_RGRORATE_m_per_y, # vertical root growth rate in m per y (accessing new soil layers)
        p_THICK,      final_Rootden_profile)

    NLAYER = length(p_THICK)
    # Hammel and Kennel 2001:

    # 0) preprocess final equilibrium state: (final_Rootden_profile)
    # final_Rootden_profile: relative values of final root density per unit volume
    # make root density continuous between first and last layer
    i1 = findfirst(final_Rootden_profile .> 1.e-6) # first layer where final_Rootden_profile[i] is >=1.e-6
    i2 = findlast(final_Rootden_profile .> 1.e-6)  # last layer (in 1:NLAYER) where final_Rootden_profile[i] is >=1.e-6
    if !isnothing(i1) && isnothing(i2)
        i2 = NLAYER
    end
    for i = i1:i2 # set final_Rootden_profile to minimum of 1.e-6 within the root zone (between i1 and i2)
        final_Rootden_profile[i] = max( final_Rootden_profile[i], 1.01e-6)
    end

    # 1) Time dependent root density parameters is a vector quantity that is dependent on time (i.e. 2D array):
    p_RELDEN_2Darray = fill(NaN, length(AGE_at_timepoints), NLAYER)

    # compute a time dependent profile if growth period RGROPER_y is not zero
    if p_RGROPER_y > 0 #zero(RGROPER_yrs)

        # 1a) Define age at which root growth arrives in a specific depth and lateral root growth commences
        # z_center     = soil_discretization[!,"Upper_m"] - soil_discr["THICK"]/1000/2 # soil depth [m] (midpoint, i.e. average depth of layer)
        z_center_m     = -(cumsum(p_THICK) - p_THICK/2)/1000
        depfirst_m  = z_center_m[1] - p_THICK[1] / 1000.
        # Here z_center_m and depfirst_m are negative values in the soil. (z axis upward)

        # 1b) compute arrival times tini
        vᵣ = p_RGRORATE_m_per_y # growth rate (m/y)
        zᵣ₀ = p_INITRDEP # initial depth of root
        # zᵣ =  # maximal depth of root
        # root_arrival_time = (zᵣ - zᵣ₀)/vᵣ
        # tᵣ =  # time period of lateral growth (after first vertical arrival)
        # rₗ₀ = #
        # # t_ini(z) : arrival time of roots at depth z, then starts lateral growth

        t_ini = fill(NaN, NLAYER)
        for i = 1:NLAYER
            t_ini[i] = 1.e+20               # start out with "never"
            z = depfirst_m - z_center_m[i]
            roots_present = final_Rootden_profile[i] >= 1.e-6
            if roots_present
                if z <= zᵣ₀
                    t_ini[i] = 0.
                elseif z > zᵣ₀
                    if vᵣ > 0
                        t_ini[i] = (z - zᵣ₀)/vᵣ
                    end
                end
            end
        end

        # 1c) lateral growth per layer starts at t_ini and ends at t_ini + tᵣ
        for i in eachindex(timepoints)
            p_RELDEN_2Darray[i,:] =
                HammelKennel_lateral_rootgrowth(;
                    final_Rootden_profile = final_Rootden_profile,
                    INITRDEP_m        = p_INITRDEP,
                    INITRLEN_m_per_m2 = p_INITRLEN,
                    t_y               = AGE_at_timepoints[i],
                    tstart_y          = t_ini,
                    RGROPER_yrs       = p_RGROPER_y)
        end
    else
        # 1a) if no growth period provieded, assume root distribution is already at constant final profile
        p_RELDEN_2Darray = fill(NaN, length(timepoints), NLAYER)
        for i in eachindex(timepoints)
            for j = 1:NLAYER
                p_RELDEN_2Darray[i,j] = final_Rootden_profile[j]
            end
        end
    end

    # 2) In all cases provide relative density as interpolated function in depth and time (even when constant through time)
    if NLAYER > 1
        p_RELDEN =  extrapolate(interpolate((timepoints, 1:NLAYER), p_RELDEN_2Darray,
                                            (Gridded(Constant{Next}()), NoInterp()), # 1st dimension: ..., 2nd dimension NoInterp()
                                            ), Flat()) # extrapolate flat, alternative: Throw()
        extrapolate(interpolate((timepoints, 1:NLAYER), p_RELDEN_2Darray,
                                            (Gridded(Constant{Next}()), NoInterp()), # 1st dimension: ..., 2nd dimension NoInterp()
                                            ), Flat()) # extrapolate flat, alternative: Throw()
    else
        p_RELDEN_1arg = extrapolate(scale(interpolate(p_RELDEN_2Darray[:,1], (BSpline(Constant{Next}()))), timepoints) ,Flat())
        function p_RELDEN(t, ignored)
            p_RELDEN_1arg(t)
        end
    end
    return p_RELDEN
end