Extend farquhar mechanism to spatially varying

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
 --------

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

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"

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;

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);
 

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

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)
 

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)
 

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(

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)
 

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}

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.

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.

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);
 

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);
 

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);
 

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
 
 

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