reduce allocations in canopy

NEWS.md

@@ -3,6 +3,10 @@ ClimaLand.jl Release Notes
 
 main
 --------
+Improve Performance of soil canopy model
+- PR [#666]
+Add base global soil canopy run to benchmark and experiments
+- PR [#591]
 
 v0.12.4
 --------

src/integrated/soil_canopy_model.jl

@@ -342,9 +342,6 @@ function lsm_radiant_energy_fluxes!(
     ϵ_soil = ground_model.parameters.emissivity
 
     # in W/m^2
-    PAR = p.canopy.radiative_transfer.par
-    NIR = p.canopy.radiative_transfer.nir
-
     LW_d_canopy = p.scratch1
     LW_u_soil = p.scratch2
     LW_net_canopy = p.canopy.radiative_transfer.LW_n
@@ -355,24 +352,24 @@ function lsm_radiant_energy_fluxes!(
     # in total: INC - OUT = CANOPY_ABS + (1-α_soil)*CANOPY_TRANS
     # SW out  = reflected par + reflected nir
     @. SW_out =
-        energy_per_photon_NIR * N_a * p.canopy.radiative_transfer.rnir +
-        energy_per_photon_PAR * N_a * p.canopy.radiative_transfer.rpar
+        energy_per_photon_NIR * N_a * p.canopy.radiative_transfer.nir.refl +
+        energy_per_photon_PAR * N_a * p.canopy.radiative_transfer.par.refl
 
     # net canopy
     @. SW_net_canopy =
-        energy_per_photon_NIR * N_a * p.canopy.radiative_transfer.anir +
-        energy_per_photon_PAR * N_a * p.canopy.radiative_transfer.apar
+        energy_per_photon_NIR * N_a * p.canopy.radiative_transfer.nir.abs +
+        energy_per_photon_PAR * N_a * p.canopy.radiative_transfer.par.abs
 
 
     # net soil = (1-α)*trans for par and nir
     @. R_net_soil .=
         energy_per_photon_NIR *
         N_a *
-        p.canopy.radiative_transfer.tnir *
+        p.canopy.radiative_transfer.nir.trans *
         (1 - α_soil_NIR) +
         energy_per_photon_PAR *
         N_a *
-        p.canopy.radiative_transfer.tpar *
+        p.canopy.radiative_transfer.par.trans *
         (1 - α_soil_PAR)
     ϵ_canopy = p.canopy.radiative_transfer.ϵ # this takes into account LAI/SAI
     @. LW_d_canopy = ((1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4) # double checked

src/shared_utilities/drivers.jl

@@ -23,6 +23,7 @@ export AbstractAtmosphericDrivers,
     snow_precipitation,
     vapor_pressure_deficit,
     displacement_height,
+    relative_humidity,
     make_update_drivers
 
 """
@@ -648,6 +649,33 @@ function vapor_pressure_deficit(T_air, P_air, q_air, thermo_params)
     return es - ea
 end
 
+"""
+    relative_humidity(T_air, P_air, q_air, thermo_params)
+
+Computes the vapor pressure deficit for air with temperature T_air,
+pressure P_air, and specific humidity q_air, using thermo_params,
+a Thermodynamics.jl param set.
+"""
+function relative_humidity(
+    T_air::FT,
+    P_air::FT,
+    q_air::FT,
+    thermo_params,
+) where {FT}
+    es = Thermodynamics.saturation_vapor_pressure(
+        thermo_params,
+        T_air,
+        Thermodynamics.Liquid(),
+    )
+    ea = Thermodynamics.partial_pressure_vapor(
+        thermo_params,
+        P_air,
+        Thermodynamics.PhasePartition(q_air),
+    )
+    return es / ea
+end
+
+
 """
     initialize_drivers(a::PrescribedAtmosphere{FT}, coords) where {FT}
 

src/standalone/Vegetation/Canopy.jl

@@ -399,14 +399,6 @@ function ClimaLand.make_update_aux(
 
         # Other auxiliary variables being updated:
         Ra = p.canopy.autotrophic_respiration.Ra
-        APAR = p.canopy.radiative_transfer.apar
-        PAR = p.canopy.radiative_transfer.par
-        TPAR = p.canopy.radiative_transfer.tpar
-        RPAR = p.canopy.radiative_transfer.rpar
-        ANIR = p.canopy.radiative_transfer.anir
-        NIR = p.canopy.radiative_transfer.nir
-        RNIR = p.canopy.radiative_transfer.rnir
-        TNIR = p.canopy.radiative_transfer.tnir
         β = p.canopy.hydraulics.β
         medlyn_factor = p.canopy.conductance.medlyn_term
         gs = p.canopy.conductance.gs
@@ -416,6 +408,9 @@ function ClimaLand.make_update_aux(
         ψ = p.canopy.hydraulics.ψ
         ϑ_l = Y.canopy.hydraulics.ϑ_l
         fa = p.canopy.hydraulics.fa
+        inc_par = p.canopy.radiative_transfer.inc_par
+        inc_nir = p.canopy.radiative_transfer.inc_nir
+        frac_diff = p.canopy.radiative_transfer.frac_diff
 
         # Current atmospheric conditions
         θs = p.drivers.θs
@@ -451,44 +446,37 @@ function ClimaLand.make_update_aux(
             canopy.radiative_transfer.parameters.ϵ_canopy *
             (1 - exp(-(LAI + SAI))) #from CLM 5.0, Tech note 4.20
         RT = canopy.radiative_transfer
+        compute_PAR!(inc_par, RT, canopy.radiation, p, t)
+        compute_NIR!(inc_nir, RT, canopy.radiation, p, t)
         K = extinction_coeff.(G_Function, θs)
-        PAR .= compute_PAR(RT, canopy.radiation, p, t)
-        NIR .= compute_NIR(RT, canopy.radiation, p, t)
-        e_sat =
-            Thermodynamics.saturation_vapor_pressure.(
-                thermo_params,
-                T_air,
-                Ref(Thermodynamics.Liquid()),
-            )
-        e =
-            Thermodynamics.partial_pressure_vapor.(
-                thermo_params,
-                P_air,
-                PhasePartition.(q_air),
-            )
-        rel_hum = e ./ e_sat
         DOY = Dates.dayofyear(
             canopy.atmos.ref_time + Dates.Second(floor(Int64, t)),
         )
-        frac_diff = @. diffuse_fraction(DOY, T_air, PAR + NIR, rel_hum, θs)
-        absorbance_tuple = compute_absorbances(
+        @. frac_diff = diffuse_fraction(
+            DOY,
+            T_air,
+            P_air,
+            q_air,
+            p.drivers.SW_d,
+            θs,
+            thermo_params,
+        )
+
+        compute_absorbances!(
+            p,
             RT,
-            PAR ./ (energy_per_photon_PAR * N_a),
-            NIR ./ (energy_per_photon_NIR * N_a),
+            inc_par,
+            inc_nir,
             LAI,
             K,
             ground_albedo_PAR(canopy.soil_driver, Y, p, t),
             ground_albedo_NIR(canopy.soil_driver, Y, p, t),
+            energy_per_photon_PAR,
+            energy_per_photon_NIR,
+            N_a,
             θs,
             frac_diff,
         )
-        APAR .= absorbance_tuple.par.abs
-        ANIR .= absorbance_tuple.nir.abs
-        TPAR .= absorbance_tuple.par.trans
-        TNIR .= absorbance_tuple.nir.trans
-        RPAR .= absorbance_tuple.par.refl
-        RNIR .= absorbance_tuple.nir.refl
-
 
         # update plant hydraulics aux
         hydraulics = canopy.hydraulics
@@ -561,7 +549,7 @@ function ClimaLand.make_update_aux(
             Vcmax25,
             canopy.photosynthesis,
             T_canopy,
-            APAR,
+            p.canopy.radiative_transfer.par.abs,
             β,
             medlyn_factor,
             c_co2_air,

src/standalone/Vegetation/canopy_parameterizations.jl

@@ -64,77 +64,102 @@ end
         K,
         α_soil_PAR,
         α_soil_NIR,
+        energy_per_photon_PAR,
+        energy_per_photon_NIR,
+        N_a,
         _,
         _,
     )
 
 Computes the PAR and NIR absorbances, reflectances, and tranmittances
 for a canopy in the case of the 
 Beer-Lambert model. The absorbances are a function of the radiative transfer 
-model, as well as the magnitude of incident PAR and NIR radiation in moles of 
-photons, the leaf area index, the extinction coefficient, and the
+model, as well as the magnitude of incident PAR and NIR radiation in W/m^2,
+ the leaf area index, the extinction coefficient, and the
 soil albedo in the PAR and NIR bands. Returns a
 NamedTuple of NamedTuple, of the form:
 (; par = (; refl = , trans = , abs = ),  nir = (; refl = , trans = , abs = ))
 """
-function compute_absorbances(
+function compute_absorbances!(
+    p,
     RT::BeerLambertModel{FT},
     PAR,
     NIR,
     LAI,
     K,
     α_soil_PAR,
     α_soil_NIR,
+    energy_per_photon_PAR,
+    energy_per_photon_NIR,
+    N_a,
     _...,
 ) where {FT}
     RTP = RT.parameters
-    canopy_par =
-        @. plant_absorbed_pfd(RT, PAR, RTP.α_PAR_leaf, LAI, K, α_soil_PAR)
-    canopy_nir =
-        @. plant_absorbed_pfd(RT, NIR, RTP.α_NIR_leaf, LAI, K, α_soil_NIR)
-    return (; par = canopy_par, nir = canopy_nir)
+    @. p.canopy.radiative_transfer.par = plant_absorbed_pfd(
+        RT,
+        PAR / (N_a * energy_per_photon_PAR),
+        RTP.α_PAR_leaf,
+        LAI,
+        K,
+        α_soil_PAR,
+    )
+    @. p.canopy.radiative_transfer.nir = plant_absorbed_pfd(
+        RT,
+        NIR / (N_a * energy_per_photon_NIR),
+        RTP.α_NIR_leaf,
+        LAI,
+        K,
+        α_soil_NIR,
+    )
 end
 
 """
-    compute_absorbances(
+    compute_absorbances!(p,
         RT::TwoStreamModel{FT},
         PAR,
         NIR,
         LAI,
         K,
         α_soil_PAR,
         α_soil_NIR,
+        energy_per_photon_PAR,
+        energy_per_photon_NIR,
+        N_a,
         θs,
         frac_diff,
     )
 
 Computes the PAR and NIR absorbances, reflectances, and tranmittances
 for a canopy in the case of the 
 Beer-Lambert model. The absorbances are a function of the radiative transfer 
-model, as well as the magnitude of incident PAR and NIR radiation in moles of 
-photons, the leaf area index, the extinction coefficient, and the
+model, as well as the magnitude of incident PAR and NIR radiation in W/m^2, 
+the leaf area index, the extinction coefficient, and the
 soil albedo in the PAR and NIR bands. 
 
 This model also depends on the diffuse fraction and the zenith angle.
 Returns a
 NamedTuple of NamedTuple, of the form:
 (; par = (; refl = , trans = , abs = ),  nir = (; refl = , trans = , abs = ))
 """
-function compute_absorbances(
+function compute_absorbances!(
+    p,
     RT::TwoStreamModel{FT},
     PAR,
     NIR,
     LAI,
     K,
     α_soil_PAR,
     α_soil_NIR,
+    energy_per_photon_PAR,
+    energy_per_photon_NIR,
+    N_a,
     θs,
     frac_diff,
 ) where {FT}
     RTP = RT.parameters
-    canopy_par = @. plant_absorbed_pfd(
+    @. p.canopy.radiative_transfer.par = plant_absorbed_pfd(
         RT,
-        PAR,
+        PAR / (energy_per_photon_PAR * N_a),
         RTP.α_PAR_leaf,
         RTP.τ_PAR_leaf,
         LAI,
@@ -143,9 +168,9 @@ function compute_absorbances(
         α_soil_PAR,
         frac_diff,
     )
-    canopy_nir = @. plant_absorbed_pfd(
+    @. p.canopy.radiative_transfer.nir = plant_absorbed_pfd(
         RT,
-        NIR,
+        NIR / (energy_per_photon_NIR * N_a),
         RTP.α_NIR_leaf,
         RTP.τ_NIR_leaf,
         LAI,
@@ -154,7 +179,6 @@ function compute_absorbances(
         α_soil_NIR,
         frac_diff,
     )
-    return (; par = canopy_par, nir = canopy_nir)
 end
 
 """
@@ -382,11 +406,11 @@ function plant_absorbed_pfd(
 end
 
 """
-    diffuse_fraction(td::FT, T::FT, SW_IN::FT, RH::FT, θs::FT) where {FT}
+    diffuse_fraction(td::FT, T::FT, P, q, SW_IN::FT, θs::FT, thermo_params) where {FT}
 
 Computes the fraction of diffuse radiation (`diff_frac`) as a function
 of the solar zenith angle (`θs`), the total surface incident shortwave radiation (`SW_IN`),
-the air temperature (`T`), the relative humidity (`RH`), and the day of the year
+the air temperature (`T`), air pressure (`P`), specific humidity (`q`), and the day of the year
 (`td`). 
 
 See Appendix A of Braghiere, "Evaluation of turbulent fluxes of CO2, sensible heat,
@@ -403,7 +427,16 @@ can become large negative numbers.
 
 Because of that, we cap the returned value to lie within [0,1].
 """
-function diffuse_fraction(td, T::FT, SW_IN::FT, RH::FT, θs::FT) where {FT}
+function diffuse_fraction(
+    td,
+    T::FT,
+    P::FT,
+    q::FT,
+    SW_IN::FT,
+    θs::FT,
+    thermo_params,
+) where {FT}
+    RH = ClimaLand.relative_humidity(T, P, q, thermo_params)
     k₀ = FT(1370 * (1 + 0.033 * cos(2π * td / 365))) * cos(θs)
     kₜ = SW_IN / k₀
     if kₜ ≤ 0.3