From 7eb5cb7a1bc3d821f6e4720c885c027416dc7258 Mon Sep 17 00:00:00 2001
From: kmdeck <kdeck@caltech.edu>
Date: Wed, 21 Aug 2024 10:00:13 -0700
Subject: [PATCH] review comments

---
 .../standalone/Vegetation/timestep_test.jl    | 26 +++++++++----------
 src/standalone/Vegetation/canopy_energy.jl    |  4 +--
 src/standalone/Vegetation/component_models.jl | 10 ++++---
 3 files changed, 21 insertions(+), 19 deletions(-)

diff --git a/experiments/standalone/Vegetation/timestep_test.jl b/experiments/standalone/Vegetation/timestep_test.jl
index 255348abf0..aebd9529a9 100644
--- a/experiments/standalone/Vegetation/timestep_test.jl
+++ b/experiments/standalone/Vegetation/timestep_test.jl
@@ -18,7 +18,7 @@ using ClimaLand.Canopy
 using ClimaLand.Canopy.PlantHydraulics
 import ClimaLand
 import ClimaLand.Parameters as LP
-const FT = Float32;
+const FT = Float64;
 earth_param_set = LP.LandParameters(FT);
 f_root_to_shoot = FT(3.5)
 plant_ν = FT(2.46e-4) # kg/m^2
@@ -94,13 +94,12 @@ AR_model = AutotrophicRespirationModel{FT}(AR_params);
 
 f_root_to_shoot = FT(3.5)
 SAI = FT(1.0)
-RAI = FT(3 * f_root_to_shoot)
+RAI = FT(3f_root_to_shoot)
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
-
-function root_distribution(z::T) where {T}
-    rooting_depth = T(1.0)
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth))
-end;
+rooting_depth = FT(1.0)
+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
 
 K_sat_plant = FT(1.8e-6)
 ψ63 = FT(-4 / 0.0098)
@@ -132,9 +131,9 @@ plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(;
     compartment_surfaces = compartment_surfaces,
     compartment_midpoints = compartment_midpoints,
 );
-
+ac_canopy = FT(1e3)
 energy_model = ClimaLand.Canopy.BigLeafEnergyModel{FT}(
-    BigLeafEnergyParameters{FT}(FT(1e3)),
+    BigLeafEnergyParameters{FT}(ac_canopy),
 )
 
 canopy = ClimaLand.Canopy.CanopyModel{FT}(;
@@ -175,10 +174,10 @@ Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1])
 Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2])
 
 
-
-t0 = 150 * 3600 * 24.0
+seconds_per_day = 3600 * 24.0
+t0 = 150seconds_per_day
 N_days = 100
-tf = t0 + 3600 * 24 * N_days
+tf = t0 + N_days * seconds_per_day
 evaluate!(Y.canopy.energy.T, atmos.T, t0)
 set_initial_cache! = make_set_initial_cache(canopy)
 set_initial_cache!(p, Y, t0);
@@ -216,8 +215,9 @@ ref_T = [parent(ref_sol.u[k].canopy.energy.T)[1] for k in 1:length(ref_sol.t)]
 mean_err = []
 p95_err = []
 p99_err = []
-dts = [225.0, 450.0, 900.0, 1800.0, 3600.0]
+dts = [12.0, 24.0, 48.0, 100.0, 225.0, 450.0, 900.0, 1800.0, 3600.0]
 for dt in dts
+    @info dt
     saveat = Array(t0:(3 * 3600):tf)
     evaluate!(Y.canopy.energy.T, atmos.T, t0)
     updateat = Array(t0:(3600 * 3):tf)
diff --git a/src/standalone/Vegetation/canopy_energy.jl b/src/standalone/Vegetation/canopy_energy.jl
index 04be57ebe4..96d4016941 100644
--- a/src/standalone/Vegetation/canopy_energy.jl
+++ b/src/standalone/Vegetation/canopy_energy.jl
@@ -188,8 +188,8 @@ function ClimaLand.make_compute_jacobian(
         area_index = p.canopy.hydraulics.area_index
         ac_canopy = model.parameters.ac_canopy
         earth_param_set = canopy.parameters.earth_param_set
-        _T_freeze = FT(LP.T_freeze(earth_param_set))
-        _σ = FT(LP.Stefan(earth_param_set))
+        _T_freeze = LP.T_freeze(earth_param_set)
+        _σ = LP.Stefan(earth_param_set)
         @. ∂LW_n∂Tc = -2 * 4 * _σ * ϵ_c * Y.canopy.energy.T^3 # ≈ ϵ_ground = 1
         @. ∂qc∂Tc = partial_q_sat_partial_T_liq(
             p.drivers.P,
diff --git a/src/standalone/Vegetation/component_models.jl b/src/standalone/Vegetation/component_models.jl
index f936e83916..ffda43199c 100644
--- a/src/standalone/Vegetation/component_models.jl
+++ b/src/standalone/Vegetation/component_models.jl
@@ -129,8 +129,9 @@ function ClimaLand.make_compute_exp_tendency(
 )
     function compute_exp_tendency!(dY, Y, p, t)
         vars = prognostic_vars(component)
-        if vars != ()
-            getproperty(getproperty(dY, name(canopy)), name(component)) .= 0
+        dY_canopy = getproperty(dY, name(canopy))
+        if !isempty(vars)
+            getproperty(dY_canopy, name(component)) .= 0
         end
     end
     return compute_exp_tendency!
@@ -152,8 +153,9 @@ function ClimaLand.make_compute_imp_tendency(
 )
     function compute_imp_tendency!(dY, Y, p, t)
         vars = prognostic_vars(component)
-        if vars != ()
-            getproperty(getproperty(dY, name(canopy)), name(component)) .= 0
+        dY_canopy = getproperty(dY, name(canopy))
+        if !isempty(vars)
+            getproperty(dY_canopy, name(component)) .= 0
         end
 
     end