fix: Edge cases of LAI=0 and IDEPTH_cm=0

Author: Fabian Bernhard

Project.toml

@@ -1,7 +1,7 @@
 name = "LWFBrook90"
 uuid = "100deb6c-4418-4b77-ad5c-a1fadb580743"
 authors = ["Fabian Bernhard"]
-version = "0.9.1"
+version = "0.9.2"
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
@@ -11,7 +11,6 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -35,7 +34,6 @@ ComponentArrays = "0.13"
 DataFrames = "0.22, 1.1"
 DataFramesMeta = "0.6, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12, 0.13"
 DiffEqCallbacks = "2"
-DifferentialEquations = "7"
 Interpolations = "0.13, 0.14"
 Literate = "2"
 Measures = "0.3"

src/LWFBrook90.jl

@@ -181,6 +181,7 @@ include("../examples/func_run_example.jl") # defines RelativeDaysFloat2DateTime
                 soil_output_depths_m::Vector = zeros(Float64, 0),
                 kwargs...)
 or
+
     remakeSPAC(parametrizedSPAC::SPAC;
                 requested_tspan = nothing,
                 soil_output_depths_m::Vector = zeros(Float64, 0),

src/func_discretize_soil_domain.jl

@@ -225,8 +225,8 @@ function refine_soil_discretization(
     # 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 Bottom_m is below IDEPTH_m
-    QDEPTH_idx = sum(is_SRFL_layer_BOOLEAN)         # lowest node where Bottom_m is below QDEPTH_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"])
@@ -607,7 +607,7 @@ function HammelKennel_transient_root_density(;
         end
     end
 
-    # 2) In all cases provide relative density as interpolated funciton in depth and time (even when constant through time)
+    # 2) In all cases provide relative density as interpolated function in depth and time (even when constant through time)
     p_fT_RELDEN =  extrapolate(interpolate((timepoints, 1:NLAYER), p_RELDEN_2Darray,
                                         (Gridded(Constant{Next}()), NoInterp()), # 1st dimension: ..., 2nd dimension NoInterp()
                                         ), Flat()) # extrapolate flat, alternative: Throw()

src/func_read_inputData.jl

@@ -123,7 +123,7 @@ function loadSPAC(folder::String, prefix::String;
             "LAI"    in names(input_meteoveg) && "SAI"    in names(input_meteoveg))
             _to_use_canopy_evolution = input_meteoveg[:, [:days, :DENSEF, :HEIGHT, :LAI, :SAI]] # store DataFrame in SPAC
             # Assert validity of vegetation values
-            @assert all(_to_use_canopy_evolution.LAI .> 0) "Error in vegetation parameters: LAI must be above 0%."
+            @assert all(_to_use_canopy_evolution.LAI .>= 0) "Error in vegetation parameters: LAI must be above 0%."
         else
             error("""
                 Input_meteoveg is expected to contain one or multiple of the columns: :DENSEF, :HEIGHT, :LAI, or :SAI.
@@ -145,8 +145,8 @@ function loadSPAC(folder::String, prefix::String;
 
         _to_use_canopy_evolution = canopy_evolution # store parameter arguments in SPAC
         # Assert validity of vegetation values
-        @assert all(_to_use_canopy_evolution.LAI.LAI_perc_BtoC > 0) "Error in vegetation parameters: LAI must be above 0%."
-        @assert all(_to_use_canopy_evolution.LAI.LAI_perc_CtoB > 0) "Error in vegetation parameters: LAI must be above 0%."
+        @assert all(_to_use_canopy_evolution.LAI.LAI_perc_BtoC >= 0) "Error in vegetation parameters: LAI must be above 0%."
+        @assert all(_to_use_canopy_evolution.LAI.LAI_perc_CtoB >= 0) "Error in vegetation parameters: LAI must be above 0%."
     end
 
     # Assert validity of vegetation values
@@ -850,10 +850,9 @@ function read_path_soil_horizons(path_soil_horizons)
     ClappHornberger_expected_column_names =
         ["HorizonNr","Upper_m","Lower_m","thsat_volFrac","thetaf_volFrac","psif_kPa","bexp_","kf_mmDay","wtinf_","gravel_volFrac"]
 
-
-    if String.(propertynames(File(path_soil_horizons))) == MualVanGen_expected_column_names
+    if strip.(String.(propertynames(File(path_soil_horizons)))) == MualVanGen_expected_column_names
         FLAG_MualVanGen = 1
-    elseif String.(propertynames(File(path_soil_horizons))) == ClappHornberger_expected_column_names
+    elseif strip.(String.(propertynames(File(path_soil_horizons)))) == ClappHornberger_expected_column_names
         FLAG_MualVanGen = 0
     else
         error("""

src/module_WAT.jl

@@ -58,7 +58,7 @@ Compute fraction of infiltration to each soil layer.
 """
 function INFPAR(p_INFEXP, p_ILAYER, p_soil)
     p_INFRAC = fill(NaN, p_soil.NLAYER) # fraction of infiltration to each layer
-    if p_INFEXP <= 0
+    if (p_INFEXP <= 0 || p_ILAYER <= 0)
         p_INFRAC[1]         = 1
         p_INFRAC[2:p_soil.NLAYER] .= 0
     else