canopy cleanup (#1181)

Author: kmdeck

src/diagnostics/default_diagnostics.jl

@@ -127,7 +127,6 @@ function default_diagnostics(
             "gpp",
             "an",
             "rd",
-            "vcmax25",
             "nir",
             "anir",
             "rnir",
@@ -295,7 +294,6 @@ function default_diagnostics(
             "gpp",
             "an",
             "rd",
-            "vcmax25",
             "nir",
             "anir",
             "rnir",

src/diagnostics/land_compute_methods.jl

@@ -211,7 +211,6 @@ end
 } p.canopy.photosynthesis.GPP
 @diagnostic_compute "photosynthesis_net_leaf" Union{SoilCanopyModel, LandModel} p.canopy.photosynthesis.An
 @diagnostic_compute "respiration_leaf" Union{SoilCanopyModel, LandModel} p.canopy.photosynthesis.Rd
-@diagnostic_compute "vcmax25" Union{SoilCanopyModel, LandModel} p.canopy.photosynthesis.Vcmax25
 
 # Canopy - Radiative Transfer
 @diagnostic_compute "near_infrared_radiation_down" Union{

src/integrated/land.jl

@@ -392,8 +392,7 @@ function lsm_radiant_energy_fluxes!(
     LW_d = p.drivers.LW_d
     SW_d = p.drivers.SW_d
 
-    T_canopy =
-        ClimaLand.Canopy.canopy_temperature(canopy.energy, canopy, Y, p, t)
+    T_canopy = ClimaLand.Canopy.canopy_temperature(canopy.energy, canopy, Y, p)
 
     α_soil_PAR = p.soil.PAR_albedo
     α_soil_NIR = p.soil.NIR_albedo

src/integrated/soil_canopy_model.jl

@@ -299,8 +299,7 @@ function lsm_radiant_energy_fluxes!(
     LW_d = p.drivers.LW_d
     SW_d = p.drivers.SW_d
 
-    T_canopy =
-        ClimaLand.Canopy.canopy_temperature(canopy.energy, canopy, Y, p, t)
+    T_canopy = ClimaLand.Canopy.canopy_temperature(canopy.energy, canopy, Y, p)
 
     α_soil_PAR = p.soil.PAR_albedo
     α_soil_NIR = p.soil.NIR_albedo

src/standalone/Vegetation/Canopy.jl

@@ -437,45 +437,22 @@ function ClimaLand.make_update_aux(
         # the current time, as they depend on prescribed fields.
         set_canopy_prescribed_field!(canopy.hydraulics, p, t)
 
-        # Other auxiliary variables being updated:
-        Ra = p.canopy.autotrophic_respiration.Ra
+        # Shortcut names
         An = p.canopy.photosynthesis.An
-        GPP = p.canopy.photosynthesis.GPP
-        Rd = p.canopy.photosynthesis.Rd
         ψ = p.canopy.hydraulics.ψ
         ϑ_l = Y.canopy.hydraulics.ϑ_l
         fa = p.canopy.hydraulics.fa
         par_d = p.canopy.radiative_transfer.par_d
         nir_d = p.canopy.radiative_transfer.nir_d
-
-        bc = canopy.boundary_conditions
-        # Current atmospheric conditions
         cosθs = p.drivers.cosθs
-        c_co2_air = p.drivers.c_co2
-        P_air = p.drivers.P
-        T_air = p.drivers.T
-        q_air = p.drivers.q
-
-        # unpack parameters
-        earth_param_set = canopy.parameters.earth_param_set
-        c = LP.light_speed(earth_param_set)
-        planck_h = LP.planck_constant(earth_param_set)
-        N_a = LP.avogadro_constant(earth_param_set)
-        grav = LP.grav(earth_param_set)
-        ρ_l = LP.ρ_cloud_liq(earth_param_set)
-        R = LP.gas_constant(earth_param_set)
-        T_freeze = LP.T_freeze(earth_param_set)
-        thermo_params = earth_param_set.thermo_params
-        (; G_Function, Ω, λ_γ_PAR) = canopy.radiative_transfer.parameters
-        energy_per_mole_photon_par = planck_h * c / λ_γ_PAR * N_a
-        (; g1, g0, Drel) = canopy.conductance.parameters
         area_index = p.canopy.hydraulics.area_index
         LAI = area_index.leaf
         SAI = area_index.stem
-        RAI = area_index.root
-        (; sc, pc) = canopy.photosynthesis.parameters
+
+        bc = canopy.boundary_conditions
 
         # update radiative transfer
+        (; G_Function, Ω, λ_γ_PAR) = canopy.radiative_transfer.parameters
         @. p.canopy.radiative_transfer.ϵ =
             canopy.radiative_transfer.parameters.ϵ_canopy *
             (1 - exp(-(LAI + SAI))) #from CLM 5.0, Tech note 4.20
@@ -554,71 +531,23 @@ function ClimaLand.make_update_aux(
                     ),
                 ) * PlantHydraulics.harmonic_mean(areaip1, areai)
         end
-        # We update the fa[n_stem+n_leaf] element once we have computed transpiration, below
-        # update photosynthesis and conductance terms
-        # This should still use T_air, P_air, q_air
-        medlyn_factor =
-            @. lazy(medlyn_term(g1, T_air, P_air, q_air, thermo_params))
-        # Anywhere we use an Arrhenius factor, use T_canopy instead T_air
-        T_canopy = canopy_temperature(canopy.energy, canopy, Y, p, t)
-
-        # update moisture stress
-        i_end = n_stem + n_leaf
-        β = @. lazy(moisture_stress(ψ.:($$i_end) * ρ_l * grav, sc, pc))
-
-        # Update Rd, An, Vcmax25 (if applicable to model) in place
-        Vcmax25 = p.canopy.photosynthesis.Vcmax25
-        update_photosynthesis!(
-            Rd,
-            An,
-            Vcmax25,
-            canopy.photosynthesis,
-            T_canopy,
-            p.canopy.radiative_transfer.par.abs,
-            β,
-            medlyn_factor,
-            c_co2_air,
-            R,
-            energy_per_mole_photon_par,
-            par_d,
-        )
+        # We update the fa[n_stem+n_leaf] element once we have computed transpiration
+
+        # Update Rd, An, Vcmax25 (if applicable to model) in place, GPP
+        update_photosynthesis!(p, Y, canopy.photosynthesis, canopy)
+
         # update SIF
-        SIF = p.canopy.sif.SIF
-        update_SIF!(
-            SIF,
-            canopy.sif,
-            p.canopy.radiative_transfer.par.abs,
-            T_canopy,
-            Vcmax25,
-            R,
-            T_freeze,
-            canopy.photosynthesis.parameters,
-            energy_per_mole_photon_par,
-            par_d,
-        )
-        @. GPP = compute_GPP(An, extinction_coeff(G_Function, cosθs), LAI, Ω)
-        @. p.canopy.conductance.r_stomata_canopy =
-            1 / upscale_leaf_conductance(
-                medlyn_conductance(g0, Drel, medlyn_factor, An, c_co2_air), #conductance, leaf level
-                LAI,
-                T_air,
-                R,
-                P_air,
-            )
+        update_SIF!(p, Y, canopy.sif, canopy)
+
+        # update stomatal conductance
+        update_canopy_conductance!(p, Y, canopy.conductance, canopy)
+
         # update autotrophic respiration
-        h_canopy = hydraulics.compartment_surfaces[end]
-        @. Ra = compute_autrophic_respiration(
+        update_autotrophic_respiration!(
+            p,
+            Y,
             canopy.autotrophic_respiration,
-            Vcmax25,
-            LAI,
-            SAI,
-            RAI,
-            extinction_coeff(G_Function, cosθs),
-            Ω,
-            An,
-            Rd,
-            β,
-            h_canopy,
+            canopy,
         )
     end
     return update_aux!

src/standalone/Vegetation/autotrophic_respiration.jl

@@ -57,6 +57,59 @@ ClimaLand.auxiliary_types(model::AutotrophicRespirationModel{FT}) where {FT} =
     (FT,)
 ClimaLand.auxiliary_domain_names(::AutotrophicRespirationModel) = (:surface,)
 
+
+
+"""
+    update_autotrophic_respiration!(p, Y, autotrophic_respiration::AutotrophicRespirationModel, canopy)
+
+Computes the autotrophic respiration rate  (mol co2 m^-2 s^-1) as the sum of the plant maintenance
+and growth respirations, according to the JULES model.
+
+Clark, D. B., et al. "The Joint UK Land Environment Simulator (JULES), model 
+description–Part 2: carbon fluxes and vegetation dynamics." Geoscientific Model Development 4.3 (2011): 701-722.
+"""
+function update_autotrophic_respiration!(
+    p,
+    Y,
+    autotrophic_respiration::AutotrophicRespirationModel,
+    canopy,
+)
+    hydraulics = canopy.hydraulics
+    n_stem = hydraulics.n_stem
+    n_leaf = hydraulics.n_leaf
+    h_canopy = hydraulics.compartment_surfaces[end]
+    i_end = n_stem + n_leaf
+    ψ = p.canopy.hydraulics.ψ
+    area_index = p.canopy.hydraulics.area_index
+    LAI = area_index.leaf
+    SAI = area_index.stem
+    RAI = area_index.root
+    earth_param_set = canopy.parameters.earth_param_set
+    grav = LP.grav(earth_param_set)
+    ρ_l = LP.ρ_cloud_liq(earth_param_set)
+    (; sc, pc) = canopy.photosynthesis.parameters
+    (; G_Function, Ω) = canopy.radiative_transfer.parameters
+    cosθs = p.drivers.cosθs
+    An = p.canopy.photosynthesis.An
+    Rd = p.canopy.photosynthesis.Rd
+
+    β = @. lazy(moisture_stress(ψ.:($$i_end) * ρ_l * grav, sc, pc))
+    Vcmax25 = get_Vcmax25(p, canopy.photosynthesis)
+    @. p.canopy.autotrophic_respiration.Ra = compute_autrophic_respiration(
+        autotrophic_respiration,
+        Vcmax25,
+        LAI,
+        SAI,
+        RAI,
+        extinction_coeff(G_Function, cosθs),
+        Ω,
+        An,
+        Rd,
+        β,
+        h_canopy,
+    )
+end
+
 """
     compute_autrophic_respiration(model::AutotrophicRespirationModel,
                                   Vcmax25,

src/standalone/Vegetation/canopy_boundary_fluxes.jl

@@ -121,7 +121,7 @@ A helper function which returns the temperature for the canopy
 model.
 """
 function ClimaLand.surface_temperature(model::CanopyModel, Y, p, t)
-    return canopy_temperature(model.energy, model, Y, p, t)
+    return canopy_temperature(model.energy, model, Y, p)
 end
 
 """

src/standalone/Vegetation/canopy_energy.jl

@@ -36,13 +36,13 @@ No equation for the energy of the canopy is solved.
 struct PrescribedCanopyTempModel{FT} <: AbstractCanopyEnergyModel{FT} end
 
 """
-    canopy_temperature(model::PrescribedCanopyTempModel, canopy, Y, p, t)
+    canopy_temperature(model::PrescribedCanopyTempModel, canopy, Y, p)
 
 Returns the canopy temperature under the `PrescribedCanopyTemp` model,
 where the canopy temperature is assumed to be the same as the atmosphere
 temperature.
 """
-function canopy_temperature(model::PrescribedCanopyTempModel, canopy, Y, p, t)
+function canopy_temperature(model::PrescribedCanopyTempModel, canopy, Y, p)
     p.drivers.T
 end
 
@@ -86,13 +86,12 @@ ClimaLand.prognostic_types(model::BigLeafEnergyModel{FT}) where {FT} = (FT,)
 ClimaLand.prognostic_domain_names(model::BigLeafEnergyModel) = (:surface,)
 
 """
-    canopy_temperature(model::BigLeafEnergyModel, canopy, Y, p, t)
+    canopy_temperature(model::BigLeafEnergyModel, canopy, Y, p)
 
 Returns the canopy temperature under the `BigLeafEnergyModel` model,
 where the canopy temperature is modeled prognostically.
 """
-canopy_temperature(model::BigLeafEnergyModel, canopy, Y, p, t) =
-    Y.canopy.energy.T
+canopy_temperature(model::BigLeafEnergyModel, canopy, Y, p) = Y.canopy.energy.T
 
 function make_compute_imp_tendency(
     model::BigLeafEnergyModel{FT},
@@ -110,19 +109,17 @@ function make_compute_imp_tendency(
         # area index on the LHF, as ac_canopy [J/m^2/K]
         # is per unit area plant.
 
+        # d(energy.T) = - net_energy_flux / specific_heat_capacity
         # To prevent dividing by zero, change AI" to
         # "max(AI, eps(FT))"
-        c_per_ground_area =
-            @. ac_canopy * max(area_index.leaf + area_index.stem, eps(FT))
 
-        # d(energy.T) = - net_energy_flux / specific_heat_capacity
         @. dY.canopy.energy.T =
             -(
                 -p.canopy.radiative_transfer.LW_n -
                 p.canopy.radiative_transfer.SW_n +
                 p.canopy.turbulent_fluxes.shf +
                 p.canopy.turbulent_fluxes.lhf - p.canopy.energy.fa_energy_roots
-            ) / c_per_ground_area
+            ) / (ac_canopy * max(area_index.leaf + area_index.stem, eps(FT)))
     end
     return compute_imp_tendency!
 end