Extend photosynthesis mechanism to support fields

This change is applied to both OptimalityFarquharParameters and FarquharParameters. Previously, the mechanism parameter was held as an empty struct of type ::C3 or ::C4 to allow for dispatch. Because it is not possible to create a ClimaCore field of different structs, the mechanism param is changed to a is_c3 param that holds 1s and 0s. This choice was made because reading a NetCDF file requires regridding, which must be performed over Numbers. If the input data file must store the mechanism as a float, then it is easiest to also internally represent the mechanism as a float. The C3/C4 structs are now removed. The functions that previously used C3/C4 structs for dispatching are converted to two seperate functions, which are called based on the is_c3 value. The constructor checks that the mechanism value of field is between 0 and 1, and throws an error otherwise. Regridding a field of 0s and 1s may result in values between 0 and 1, so the constuctor rounds them to the nearest integer. This implementation allows for future support of mixed mechanisms in an area represented by proportion c3 without many changes. All calls to the FarquharParameters constructor are changed to work with the changes to mechanism representation. Tests are also changed to work with new constructor and representation of mechanism. The tests in canopy_model are changed to run twice using a for loop. The first run has all the parameters as floats, and the second run has all canopy parameters that can be fields as fields. In this case, the photosynthesis mechanism is tested as a field that alternates between 1.0 and 0.0, which represents alternating between C3 and C4. The tests were previously all run over a Point domain. The tests that use FarquharParameters are changed to run on a spherical surface. Additionally, a test for vcmax as a field is added. Land.jl and land_region.jl are changed to use the mechanism map from the clm data.
CliMA · Sep 17, 2024 · 3e9a73a · 3e9a73a
1 parent b47ef0f
commit 3e9a73a
Show file tree

Hide file tree

Showing 28 changed files with 244 additions and 159 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -4,6 +4,8 @@ ClimaLand.jl Release Notes
 main
 --------
  - Add regional simulation example PR [#757](https://github.com/CliMA/ClimaLand.jl/pull/757)
+ - ![][badge-💥breaking] Extend photosynthesis mechanism parameter to support fields.
+  PR[#774](https://github.com/CliMA/ClimaLand.jl/pull/774)
 
 v0.14.3
 --------

diff --git a/docs/src/standalone/apps/leaf_an.jl b/docs/src/standalone/apps/leaf_an.jl
@@ -38,10 +38,11 @@ function ParamViz.parameterisation(
     Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R)
     medlynt = 1 + g1 / sqrt(VPD)
     ci = intercellular_co2(ca, Γstar, medlynt)
-    Aj = light_assimilation(Canopy.C3(), J, ci, Γstar)
+    is_c3 = 1.0
+    Aj = light_assimilation(is_c3, J, ci, Γstar)
     Kc = MM_Kc(Kc25, ΔHkc, T, To, R)
     Ko = MM_Ko(Ko25, ΔHko, T, To, R)
-    Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi)
+    Ac = rubisco_assimilation(is_c3, Vcmax, ci, Γstar, Kc, Ko, oi)
     Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R)
     An = net_photosynthesis(Ac, Aj, Rd, β)
     return An
@@ -83,7 +84,7 @@ function An_app_f()
             "Ω", # Clumping index
             "Γstar25", # co2 compensation at 25c
             "ΔHΓstar", # a constant energy of activation
-            "To", # a standard temperature 
+            "To", # a standard temperature
             "R", # the universal gas constant
             "Vcmax25", # the maximum rate of carboxylation of Rubisco
             "ΔHJmax", # a constant

diff --git a/docs/src/standalone/apps/leaf_an_ci.jl b/docs/src/standalone/apps/leaf_an_ci.jl
@@ -36,8 +36,9 @@ function ParamViz.parameterisation(
     Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R)
     Kc = MM_Kc(Kc25, ΔHkc, T, To, R)
     Ko = MM_Ko(Ko25, ΔHko, T, To, R)
-    Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi)
-    Aj = light_assimilation(Canopy.C3(), J, ci, Γstar)
+    is_c3 = 1.0
+    Ac = rubisco_assimilation(is_c3, Vcmax, ci, Γstar, Kc, Ko, oi)
+    Aj = light_assimilation(is_c3, J, ci, Γstar)
     Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R)
     An = net_photosynthesis(Ac, Aj, Rd, β)
     return An
@@ -79,7 +80,7 @@ function An_ci_app_f()
             "Ω", # Clumping index
             "Γstar25", # co2 compensation at 25c
             "ΔHΓstar", # a constant energy of activation
-            "To", # a standard temperature 
+            "To", # a standard temperature
             "R", # the universal gas constant
             "ΔHJmax", # a constant
             "θj", # an empirical "curvature parameter"

diff --git a/docs/tutorials/integrated/soil_canopy_tutorial.jl b/docs/tutorials/integrated/soil_canopy_tutorial.jl
@@ -231,8 +231,9 @@ radiative_transfer_args = (;
 
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1 = 141))
 
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5)));
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5)));
 
 K_sat_plant = FT(1.8e-8)
 RAI = (SAI + maxLAI) * f_root_to_shoot;

diff --git a/docs/tutorials/standalone/Canopy/canopy_tutorial.jl b/docs/tutorials/standalone/Canopy/canopy_tutorial.jl
@@ -98,7 +98,7 @@ land_domain = Point(; z_sfc = FT(0.0))
 #   the simulation with atmospheric and radiative flux models. We also
 # read in the observed LAI and let that vary in time in a prescribed manner.
 
-# Use the data tools for reading FLUXNET data sets 
+# Use the data tools for reading FLUXNET data sets
 include(
     joinpath(pkgdir(ClimaLand), "experiments/integrated/fluxnet/data_tools.jl"),
 );
@@ -174,8 +174,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
 # Arguments for photosynthesis model:
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 

diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
@@ -380,6 +380,14 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         regridder_kwargs = (; extrapolation_bc,),
         file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
     )
+    # photosynthesis mechanism is read as a float, where 1.0 indicates c3 and 0.0 c4
+    is_c3 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "mechanism_map.nc"),
+        "c3_dominant",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
 
     # Plant Hydraulics and general plant parameters
     SAI = FT(0.0) # m2/m2
@@ -465,13 +473,8 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     conductance_args =
         (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
     # Set up photosynthesis
-    photosynthesis_args = (;
-        parameters = Canopy.FarquharParameters(
-            FT,
-            Canopy.C3();
-            Vcmax25 = Vcmax25,
-        )
-    )
+    photosynthesis_args =
+        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
     # Set up plant hydraulics
     # Note that we clip all values of LAI below 0.05 to zero.
     # This is because we currently run into issues when LAI is

diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -184,8 +184,9 @@ radiative_transfer_args = (;
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 

diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -141,8 +141,9 @@ radiative_transfer_args = (;
 # Set up conductance
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 

diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl
@@ -222,9 +222,9 @@ radiative_transfer_args = (;
 # Set up conductance
 conductance_args = (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
-photosynthesis_args = (;
-    parameters = Canopy.FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25)
-)
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photosynthesis_args =
+    (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 # Not ideal
 LAIfunction = TimeVaryingInput(

diff --git a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
@@ -129,8 +129,9 @@ radiative_transfer_args = (;
 # Set up conductance
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 

diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl
@@ -266,8 +266,9 @@ radiative_transfer_args = (; parameters = TwoStreamParameters(FT))
 # Set up conductance
 conductance_args = (; parameters = MedlynConductanceParameters(FT; g1))
 # Set up photosynthesis
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
 photosynthesis_args =
-    (; parameters = FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25))
+    (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 function root_distribution(z::T) where {T}

diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
@@ -389,7 +389,14 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         regridder_kwargs = (; extrapolation_bc,),
         file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
     )
-
+    # photosynthesis mechanism is read as a float, where 1.0 indicates c3 and 0.0 c4
+    is_c3 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "mechanism_map.nc"),
+        "c3_dominant",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
     # Plant Hydraulics and general plant parameters
     SAI = FT(0.0) # m2/m2
     f_root_to_shoot = FT(3.5)
@@ -474,13 +481,8 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     conductance_args =
         (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
     # Set up photosynthesis
-    photosynthesis_args = (;
-        parameters = Canopy.FarquharParameters(
-            FT,
-            Canopy.C3();
-            Vcmax25 = Vcmax25,
-        )
-    )
+    photosynthesis_args =
+        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
     # Set up plant hydraulics
 
     # Note that we clip all values of LAI below 0.05 to zero.

diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl
@@ -368,8 +368,18 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
 
     # Energy Balance model
     ac_canopy = FT(2.5e4) # this will likely be 10x smaller!
-
     clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
+    # Conductance Model
+    # g1 is read in units of sqrt(kPa) and then converted to sqrt(Pa)
+    g1 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "medlynslope",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+        file_reader_kwargs = (; preprocess_func = (data) -> data * 10^(3 / 2),),
+    )
+
     # vcmax is read in units of umol CO2/m^2/s and then converted to mol CO2/m^2/s
     Vcmax25 = SpaceVaryingInput(
         joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
@@ -380,14 +390,13 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         file_reader_kwargs = (; preprocess_func = (data) -> data / 1_000_000,),
     )
 
-    # g1 is read in units of sqrt(kPa) and then converted to sqrt(Pa)
-    g1 = SpaceVaryingInput(
-        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
-        "medlynslope",
+    # c3 is read in as a float
+    is_c3 = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "mechanism_map.nc"),
+        "c3_dominant",
         surface_space;
         regridder_type,
         regridder_kwargs = (; extrapolation_bc,),
-        file_reader_kwargs = (; preprocess_func = (data) -> data * 10^(3 / 2),),
     )
 
     # Plant Hydraulics and general plant parameters
@@ -474,13 +483,8 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
     conductance_args =
         (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
     # Set up photosynthesis
-    photosynthesis_args = (;
-        parameters = Canopy.FarquharParameters(
-            FT,
-            Canopy.C3();
-            Vcmax25 = Vcmax25,
-        )
-    )
+    photosynthesis_args =
+        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
     # Set up plant hydraulics
 
     # Note that we clip all values of LAI below 0.05 to zero.

diff --git a/experiments/standalone/Vegetation/no_vegetation.jl b/experiments/standalone/Vegetation/no_vegetation.jl
@@ -82,8 +82,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 

diff --git a/experiments/standalone/Vegetation/varying_lai.jl b/experiments/standalone/Vegetation/varying_lai.jl
@@ -82,8 +82,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 

diff --git a/experiments/standalone/Vegetation/varying_lai_with_stem.jl b/experiments/standalone/Vegetation/varying_lai_with_stem.jl
@@ -82,8 +82,8 @@ cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
 
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
-
-photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
+is_c3 = FT(1) # set the photosynthesis mechanism to C3
+photo_params = FarquharParameters(FT, is_c3; Vcmax25 = FT(5e-5))
 
 photosynthesis_model = FarquharModel{FT}(photo_params);
 

diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl
@@ -107,14 +107,14 @@ end
 """
     function FarquharParameters(
         FT,
-        mechanism::AbstractPhotosynthesisMechanism;
+        is_c3::Union{FT, ClimaCore.Fields.Field};
         Vcmax25 = FT(5e-5),
         kwargs...  # For individual parameter overrides
     )
 
     function FarquharParameters(
         toml_dict::CP.AbstractTOMLDict,
-        mechanism::AbstractPhotosynthesisMechanism;
+        is_c3::Union{AbstractFloat, ClimaCore.Fields.Field};
         Vcmax25 = FT(5e-5),
         kwargs...  # For individual parameter overrides
     )
@@ -123,25 +123,25 @@ Constructors for the FarquharParameters struct. Two variants:
 1. Pass in the float-type and retrieve parameter values from the default TOML dict.
 2. Pass in a TOML dictionary to retrieve parameter values.Possible calls:
 ```julia
-ClimaLand.Canopy.FarquharParameters(Float64, ClimaLand.Canopy.C3())
+ClimaLand.Canopy.FarquharParameters(Float64, 1.0)
 # Kwarg overrides
-ClimaLand.Canopy.FarquharParameters(Float64, ClimaLand.Canopy.C3(); Vcmax25 = 99999999, pc = 444444444)
+ClimaLand.Canopy.FarquharParameters(Float64, 1.0; Vcmax25 = 99999999, pc = 444444444)
 # TOML Dictionary:
 import ClimaParams as CP
 toml_dict = CP.create_toml_dict(Float32);
-ClimaLand.Canopy.FarquharParameters(toml_dict, ClimaLand.Canopy.C3(); Vcmax25 = 99999999, pc = 444444444)
+ClimaLand.Canopy.FarquharParameters(toml_dict, 1.0f0; Vcmax25 = 99999999, pc = 444444444)
 ```
 """
 FarquharParameters(
     ::Type{FT},
-    mechanism;
+    is_c3::Union{FT, ClimaCore.Fields.Field};
     kwargs...,
 ) where {FT <: AbstractFloat} =
-    FarquharParameters(CP.create_toml_dict(FT), mechanism; kwargs...)
+    FarquharParameters(CP.create_toml_dict(FT), is_c3; kwargs...)
 
 function FarquharParameters(
     toml_dict::CP.AbstractTOMLDict,
-    mechanism;
+    is_c3::Union{AbstractFloat, ClimaCore.Fields.Field};
     Vcmax25 = 5e-5,
     kwargs...,
 )
@@ -165,11 +165,24 @@ function FarquharParameters(
     )
     parameters = CP.get_parameter_values(toml_dict, name_map, "Land")
     FT = CP.float_type(toml_dict)
-    MECH = typeof(mechanism)
+    if maximum(is_c3) > 1
+        error(
+            "is_c3 has maximum of $(maximum(is_c3)). is_c3 should be between 0 and 1",
+        )
+    end
+    if minimum(is_c3) < 0
+        error(
+            "is_c3 has minimum of $(minimum(is_c3)). is_c3 should be between 0 and 1",
+        )
+    end
+    # if is_c3 is a field, is_c3 may contain values between 0.0 and 1.0 after regridding
+    # this deals with that possibility by rounding to the closest int
+    is_c3 = round.(is_c3)
+    MECH = typeof(is_c3)
     Vcmax25 = FT.(Vcmax25)
     VC = typeof(Vcmax25)
     return FarquharParameters{FT, MECH, VC}(;
-        mechanism,
+        is_c3,
         Vcmax25,
         parameters...,
         kwargs...,
@@ -229,8 +242,8 @@ function OptimalityFarquharParameters(toml_dict; kwargs...)
 
     params = CP.get_parameter_values(toml_dict, name_map, "Land")
     FT = CP.float_type(toml_dict)
-    mechanism = ClimaLand.Canopy.C3()
-    return OptimalityFarquharParameters{FT}(; params..., kwargs..., mechanism)
+    is_c3 = FT(1)
+    return OptimalityFarquharParameters{FT, FT}(; params..., kwargs..., is_c3)
 end
 
 

diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl
@@ -139,7 +139,7 @@ function conductance_harvard(;
 end
 
 function photosynthesis_harvard(;
-    mechanism = Canopy.C3(),
+    is_c3 = FT(1),
     oi = FT(0.209),
     ϕ = FT(0.6),
     θj = FT(0.9),
@@ -176,7 +176,7 @@ function photosynthesis_harvard(;
         f,
         sc,
         pc,
-        mechanism,
+        is_c3,
     )
 end