area index in aux

still wip

Tests pass

gitignore update

Start to work on rev comments

renaming

file rename

lai consistency check

addressing Charlie's comment

revert rename of file

fixed docs

added a test

fixed test

CC 10.44

fixing tests again

docs bug

.gitignore

@@ -5,6 +5,7 @@
 *~
 *.png
 *.pdf
+*.jld2
 *.csv
 *.g
 # Files generated by invoking Julia with --code-coverage

Project.toml

@@ -25,7 +25,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 ArtifactWrappers = "0.2"
 CFTime = "0.1"
 ClimaComms = "0.3, 0.4, 0.5"
-ClimaCore = "0.10.32"
+ClimaCore = "0.10.44"
 ClimaCoreTempestRemap = "0.3.5"
 DiffEqCallbacks = "2"
 DocStringExtensions = "0.8, 0.9"

docs/tutorials/Canopy/canopy_tutorial.jl

@@ -100,14 +100,10 @@ include(
 
 # Populate the SharedCanopyParameters struct, which holds the parameters 
 # shared between all different components of the canopy model.
-
-LAI = FT(4.2)
 z0_m = FT(2)
 z0_b = FT(0.2)
 
 shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
-    LAI,
-    h_stem + h_leaf,
     z0_m,
     z0_b,
     earth_param_set,
@@ -162,13 +158,15 @@ photosynthesis_model = FarquharModel{FT}(photo_params);
 
 # Arguments for plant hydraulics model are more complicated. 
 
-# Begin by providing general plant parameters:
-
-
+# Begin by providing general plant parameters. For the area
+#indices of the canopy, we choose a `PrescribedSiteAreaIndex`,
+# which supports LAI as a function of time, with RAI and SAI as constant.
+LAI = 4.2
+LAIfunction = (t) -> eltype(t)(LAI)
 SAI = FT(0.00242)
 f_root_to_shoot = FT(3.5)
 RAI = (SAI + LAI) * f_root_to_shoot
-area_index = (root = RAI, stem = SAI, leaf = LAI)
+ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 rooting_depth = FT(1.0);
 
 # Define the root distribution function p(z):
@@ -197,7 +195,7 @@ retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a);
 S_s = FT(1e-2 * 0.0098)
 
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
-    area_index = area_index,
+    ai_parameterization = ai_parameterization,
     ν = ν,
     S_s = S_s,
     root_distribution = root_distribution,
@@ -264,12 +262,6 @@ for i in 1:2
     Y.canopy.hydraulics.ϑ_l[i] .= augmented_liquid_fraction.(ν, S_l_ini[i])
 end;
 
-# Initialize the auxiliary variables for the canopy using the initial 
-# conditions
-
-update_aux! = make_update_aux(canopy)
-update_aux!(p, Y, 0.0);
-
 # Select a time range to perform time stepping over, and a dt. Also create the 
 # saveat Array to contain the data from the model at each time step. As usual,
 # the timestep depends on the problem you are solving, the accuracy of the
@@ -279,6 +271,16 @@ t0 = FT(0)
 N_days = 365
 tf = t0 + FT(3600 * 24 * N_days)
 dt = FT(225)
+
+# Initialize the auxiliary variables for the canopy using the initial 
+# conditions and initial time.
+
+set_initial_aux_state! = make_set_initial_aux_state(canopy)
+set_initial_aux_state!(p, Y, t0);
+
+# Allocate the struct which stores the saved auxiliary state
+# and create the callback which saves it at each element in saveat.
+
 n = 16
 saveat = Array(t0:(n * dt):tf)
 

docs/tutorials/Canopy/soil_plant_hydrology_tutorial.jl

@@ -169,11 +169,10 @@ photosynthesis_args = (;
 f_root_to_shoot = FT(3.5)
 SAI = FT(0.00242)
 LAI = FT(4.2)
+LAIfunction = (t) -> eltype(t)(LAI)
 K_sat_plant = 1.8e-8
 RAI = (SAI + LAI) * f_root_to_shoot
-
-area_index = (root = RAI, stem = SAI, leaf = LAI)
-
+ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 function root_distribution(z::T; rooting_depth = FT(1.0)) where {T}
     return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
 end
@@ -191,7 +190,7 @@ plant_ν = FT(0.7)
 plant_S_s = FT(1e-2 * 0.0098)
 
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
-    area_index = area_index,
+    ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
     root_distribution = root_distribution,
@@ -224,15 +223,11 @@ canopy_component_args = (;
     hydraulics = plant_hydraulics_args,
 );
 
-# We also need to provide the shared parameter struct to the canopy as usual and
-# specify the domain for the canopy.
-
+# We also need to provide the shared parameter struct to the canopy.
 z0_m = FT(2)
 z0_b = FT(0.2)
 
 shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
-    LAI,
-    h_stem + h_leaf,
     z0_m,
     z0_b,
     earth_param_set,
@@ -304,8 +299,8 @@ end;
 # Now initialize the auxiliary variables for the combined soil and plant model.
 
 t0 = FT(0)
-update_aux! = make_update_aux(land)
-update_aux!(p, Y, t0);
+set_initial_aux_state! = make_set_initial_aux_state(land)
+set_initial_aux_state!(p, Y, t0);
 
 # Select the timestepper and solvers needed for the specific problem, and setup
 # the jacobian for the soil model. For more details on the soil jacobian setup 

experiments/integrated/ozark/hydrology_only/ozark.jl

@@ -102,14 +102,14 @@ photosynthesis_args = (;
     )
 )
 # Set up plant hydraulics
-area_index = (root = RAI, stem = SAI, leaf = LAI)
+ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
 function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
     return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
 end
 
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
-    area_index = area_index,
+    ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
     root_distribution = root_distribution,
@@ -134,8 +134,6 @@ canopy_component_args = (;
 )
 # Other info needed
 shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
-    LAI,
-    h_stem + h_leaf,
     z0_m,
     z0_b,
     earth_param_set,
@@ -179,8 +177,8 @@ for i in 1:2
         augmented_liquid_fraction.(plant_ν, S_l_ini[i])
 end
 
-update_aux! = make_update_aux(land)
-update_aux!(p, Y, t0)
+set_initial_aux_state! = make_set_initial_aux_state(land)
+set_initial_aux_state!(p, Y, t0);
 
 # Set up timestepper and jacobian for soil
 update_jacobian! = make_update_jacobian(land.soil)

experiments/integrated/ozark/ozark.jl

@@ -111,14 +111,14 @@ photosynthesis_args = (;
     )
 )
 # Set up plant hydraulics
-area_index = (root = RAI, stem = SAI, leaf = LAI)
+ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
 function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
     return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
 end
 
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
-    area_index = area_index,
+    ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
     root_distribution = root_distribution,
@@ -142,8 +142,6 @@ canopy_component_args = (;
 )
 # Other info needed
 shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
-    LAI,
-    h_stem + h_leaf,
     z0_m,
     z0_b,
     earth_param_set,
@@ -196,8 +194,8 @@ for i in 1:2
         augmented_liquid_fraction.(plant_ν, S_l_ini[i])
 end
 
-update_aux! = make_update_aux(land)
-update_aux!(p, Y, t0)
+set_initial_aux_state! = make_set_initial_aux_state(land)
+set_initial_aux_state!(p, Y, t0);
 
 # Simulation
 sv = (;

experiments/integrated/ozark/ozark_parameters.jl

@@ -60,6 +60,7 @@ To = FT(298.15)
 # Plant Hydraulics and general plant parameters
 SAI = FT(1.0) # m2/m2 or: estimated from Wang et al, FT(0.00242) ?
 LAI = FT(4.2) # m2/m2, from Wang et al.
+LAIfunction = (t) -> eltype(t)(LAI)
 f_root_to_shoot = FT(3.5)
 RAI = (SAI + LAI) * f_root_to_shoot # CLM
 K_sat_plant = 5e-9 # m/s # seems much too small?

src/ClimaLSM.jl

@@ -25,7 +25,7 @@ include("shared_utilities/sources.jl")
 include("shared_utilities/implicit_tendencies.jl")
 include("shared_utilities/implicit_functions.jl")
 include("standalone/Bucket/Bucket.jl")
-export make_interactions_update_aux, domain
+export make_interactions_update_aux
 """
      AbstractLandModel{FT} <: AbstractModel{FT}
 
@@ -51,15 +51,6 @@ abstract type AbstractLandModel{FT} <: AbstractModel{FT} end
 
 ClimaLSM.name(::AbstractLandModel) = :land
 
-"""
-    domain(model::AbstractModel)
-
-Returns a symbol indicating the model's domain, e.g. :surface or :subsurface
-"""
-domain(model::AbstractModel) = error(
-    "`domain` not implemented for $(Base.typename(typeof(model)).wrapper)",
-)
-
 function Domains.coordinates(model::AbstractLandModel)
     components = land_components(model)
     coords_list = map(components) do (component)
@@ -166,6 +157,21 @@ function make_exp_tendency(land::AbstractLandModel)
     return exp_tendency!
 end
 
+function make_set_initial_aux_state(land::AbstractLandModel)
+    interactions_update_aux! = make_interactions_update_aux(land)
+    components = land_components(land)
+    # These functions also call update_aux
+    set_initial_aux_function_list =
+        map(x -> make_set_initial_aux_state(getproperty(land, x)), components)
+    function set_initial_aux_state!(p, Y0, t0)
+        for f! in set_initial_aux_function_list
+            f!(p, Y0, t0)
+        end
+        interactions_update_aux!(p, Y0, t0) # this has to come last.
+    end
+    return set_initial_aux_state!
+end
+
 function make_update_aux(land::AbstractLandModel)
     interactions_update_aux! = make_interactions_update_aux(land)
     components = land_components(land)

src/integrated/soil_canopy_model.jl

@@ -221,11 +221,14 @@ function make_interactions_update_aux(
 ) where {FT, SM <: Soil.EnergyHydrology{FT}, RM <: Canopy.CanopyModel{FT}}
     function update_aux!(p, Y, t)
         z = ClimaCore.Fields.coordinate_field(land.soil.domain.space).z
-        (; area_index, conductivity_model) = land.canopy.hydraulics.parameters
+        (; conductivity_model) = land.canopy.hydraulics.parameters
+        area_index = p.canopy.hydraulics.area_index
         @. p.root_extraction =
             (
-                area_index[:root] +
-                area_index[land.canopy.hydraulics.compartment_labels[1]]
+                area_index.root + getproperty(
+                    area_index,
+                    land.canopy.hydraulics.compartment_labels[1],
+                )
             ) / 2 *
             PlantHydraulics.flux(
                 z,

src/integrated/soil_plant_hydrology_model.jl

@@ -134,11 +134,14 @@ function make_interactions_update_aux(
 ) where {FT, SM <: Soil.RichardsModel{FT}, RM <: Canopy.CanopyModel{FT}}
     function update_aux!(p, Y, t)
         z = ClimaCore.Fields.coordinate_field(land.soil.domain.space).z
-        (; area_index, conductivity_model) = land.canopy.hydraulics.parameters
+        (; conductivity_model) = land.canopy.hydraulics.parameters
+        area_index = p.canopy.hydraulics.area_index
         @. p.root_extraction =
             (
-                area_index[:root] +
-                area_index[land.canopy.hydraulics.compartment_labels[1]]
+                area_index.root + getproperty(
+                    area_index,
+                    land.canopy.hydraulics.compartment_labels[1],
+                )
             ) / 2 *
             PlantHydraulics.flux(
                 z,

src/shared_utilities/models.jl

@@ -15,7 +15,8 @@ export AbstractModel,
     prognostic_types,
     auxiliary_types,
     make_set_initial_aux_state,
-    name
+    name,
+    domain
 
 import .Domains: coordinates
 ## Default methods for all models - to be in a seperate module at some point.
@@ -53,6 +54,17 @@ Returns a symbol of the model component name, e.g. :soil or :vegetation.
 name(model::AbstractModel) =
     error("`name` not implemented for $(Base.typename(typeof(model)).wrapper)")
 
+"""
+    domain(model::AbstractModel)
+
+Returns a symbol indicating the model's domain, e.g. :surface or :subsurface. Only required for models that will be used as part of an LSM.
+"""
+domain(model::AbstractModel) = error(
+    "`domain` not implemented for $(Base.typename(typeof(model)).wrapper)",
+)
+
+
+
 """
    prognostic_vars(m::AbstractModel)