From eb00ebaeaa444143874eb6c0ac858e70716a79bf Mon Sep 17 00:00:00 2001
From: AlexisRenchon <a.renchon@gmail.com>
Date: Tue, 16 Jul 2024 03:57:47 -0700
Subject: [PATCH] optional arg for default diagnostics

---
 src/diagnostics/default_diagnostics.jl        |  18 +-
 src/diagnostics/soilcanopy_compute_methods.jl | 192 ++++++++++++------
 2 files changed, 147 insertions(+), 63 deletions(-)

diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl
index feadb6e2f3..9eb585bd51 100644
--- a/src/diagnostics/default_diagnostics.jl
+++ b/src/diagnostics/default_diagnostics.jl
@@ -62,7 +62,7 @@ function default_diagnostics(land_model::BucketModel, t_start; output_writer)
         "wsoil",
         "wsfc",
         "ssfc",
-    ] # TO DO: would it be helpful to return this list?
+    ]
 
     default_outputs =
         hourly_averages(bucket_diagnostics...; output_writer, t_start)
@@ -75,6 +75,7 @@ function default_diagnostics(
     t_start;
     output_writer,
     output_vars = :long,
+    average_period = :hourly,
 )
 
     define_diagnostics!(land_model)
@@ -96,7 +97,7 @@ function default_diagnostics(
             "gs",
             "mt",
             "trans",
-            "rain", # do we want?
+            "rain",
             "an",
             "gpp",
             "rd",
@@ -118,8 +119,17 @@ function default_diagnostics(
         soilcanopy_diagnostics = ["gpp", "ct", "ai", "slw", "si"]
     end
 
-    default_outputs =
-        hourly_averages(soilcanopy_diagnostics...; output_writer, t_start)
+    if average_period == :hourly
+        default_outputs =
+            hourly_averages(soilcanopy_diagnostics...; output_writer, t_start)
+    elseif average_period == :daily
+        default_outputs =
+            daily_averages(soilcanopy_diagnostics...; output_writer, t_start)
+    elseif average_period == :monthly # !! this is currently 30 days, not exact month
+        default_outputs =
+            monthly_averages(soilcanopy_diagnostics...; output_writer, t_start)
+    end
+
     return [default_outputs...]
 end
 
diff --git a/src/diagnostics/soilcanopy_compute_methods.jl b/src/diagnostics/soilcanopy_compute_methods.jl
index eef836b786..388b4fa554 100644
--- a/src/diagnostics/soilcanopy_compute_methods.jl
+++ b/src/diagnostics/soilcanopy_compute_methods.jl
@@ -1,3 +1,81 @@
+macro land_compute(name, model, compute)
+    functions = Expr[]
+    function_name = Symbol("compute_", name, "!")
+    function_body = esc(quote
+        function $function_name(out, Y, p, t, land_model::$model)
+            if isnothing(out)
+                return copy($compute)
+            else
+                out .= $compute
+            end
+        end
+    end)
+    push!(functions, function_body)
+    return quote
+        $(functions...)
+    end
+end
+
+# note: model can be e.g., BucketModel, SoilCanopyModel, or Union{BucketModel, SoilCanopyModel}
+
+@land_compute "soil_net_radiation" SoilCanopyModel p.soil.R_n
+@land_compute "soil_latent_heat_flux" SoilCanopyModel p.soil.turbulent_flux.lhf
+@land_compute "soil_aerodynamic_resistance" SoilCanopyModel p.soil.turbulent_fluxes.r_ae
+@land_compute "soil_sensible_heat_flux" SoilCanopyModel p.soil.turbulent_fluxes.shf
+@land_compute "vapor_flux" SoilCanopyModel p.soil.turbulent_fluxes.vapor_flux
+@land_compute "soil_temperature" SoilCanopyModel p.soil.T
+@land_compute "soil_water_liquid" SoilCanopyModel top_center_to_surface(
+    p.soil.θ_l,
+)
+@land_compute "infiltration" SoilCanopyModel p.soil.infiltration
+@land_compute "soilco2_diffusivity" SoilCanopyModel p.soilco2.D
+@land_compute "soilco2_source_microbe" SoilCanopyModel p.soilco2.Sm
+@land_compute "stomatal_conductance" SoilCanopyModel p.canopy.conductance.gs
+@land_compute "medlyn_term" SoilCanopyModel p.canopy.conductance.medlyn_term
+@land_compute "canopy_transpiration" SoilCanopyModel p.canopy.transpiration
+@land_compute "rainfall" SoilCanopyModel p.drivers.P_liq
+@land_compute "snowfall" SoilCanopyModel p.drivers.P_snow
+@land_compute "pressure" SoilCanopyModel p.drivers.P
+@land_compute "wind_speed" SoilCanopyModel p.drivers.u
+@land_compute "specific_humidity" SoilCanopyModel p.drivers.q
+@land_compute "air_co2" SoilCanopyModel p.drivers.c_co2
+@land_compute "radiation_shortwave_down" SoilCanopyModel p.drivers.SW_d
+@land_compute "radiation_longwave_down" SoilCanopyModel p.drivers.LW_d
+@land_compute "photosynthesis_net_leaf" SoilCanopyModel p.canopy.photosynthesis.An
+@land_compute "photosynthesis_net_canopy" SoilCanopyModel p.canopy.photosynthesis.GPP
+@land_compute "respiration_leaf" SoilCanopyModel p.canopy.photosynthesis.Rd
+@land_compute "vcmax25" SoilCanopyModel p.canopy.photosynthesis.vcmax25
+@land_compute "photosynthetically_active_radiation" SoilCanopyModel p.canopy.radiative_transfer.par
+@land_compute "photosynthetically_active_radiation_absorbed" SoilCanopyModel p.canopy.radiative_transfer.apar
+@land_compute "photosynthetically_active_radiation_reflected" SoilCanopyModel p.canopy.radiative_transfer.rpar
+@land_compute "photosynthetically_active_radiation_transmitted" SoilCanopyModel p.canopy.radiative_transfer.tpar
+@land_compute "near_infrared_radiation" SoilCanopyModel p.canopy.radiative_transfer.nir
+@land_compute "near_infrared_radiation_absorbed" SoilCanopyModel p.canopy.radiative_transfer.anir
+@land_compute "near_infrared_radiation_reflected" SoilCanopyModel p.canopy.radiative_transfer.rnir
+@land_compute "near_infrared_radiation_transmitted" SoilCanopyModel p.canopy.radiative_transfer.tnir
+@land_compute "radiation_shortwave_net" SoilCanopyModel p.canopy.radiative_transfer.SW_n
+@land_compute "radiation_longwave_net" SoilCanopyModel p.canopy.radiative_transfer.LW_n
+@land_compute "autotrophic_respiration" SoilCanopyModel p.canopy.autotrophic_respiration.Ra
+@land_compute "soilco2" SoilCanopyModel Y.soilco2.C
+@land_compute "heterotrophic_respiration" SoilCanopyModel p.soilco2.top_bc
+@land_compute "soil_hydraulic_conductivity" SoilCanopyModel p.soil.K
+@land_compute "soil_water_potential" SoilCanopyModel p.soil.ψ
+@land_compute "soil_thermal_conductivity" SoilCanopyModel p.soil.κ
+@land_compute "solar_zenith_angle" SoilCanopyModel p.drivers.θs
+@land_compute "moisture_stress_factor" SoilCanopyModel p.canopy.hydraulics.β
+@land_compute "canopy_water_potential" SoilCanopyModel p.canopy.hydraulics.ψ
+@land_compute "cross_section" SoilCanopyModel p.canopy.hydraulics
+@land_compute "cross_section_roots" SoilCanopyModel p.canopy.hydraulics.fa_roots
+@land_compute "area_index" SoilCanopyModel p.canopy.hydraulics.area_index.leaf
+@land_compute "canopy_latent_heat_flux" SoilCanopyModel p.canopy.energy.lhf
+@land_compute "canopy_sensible_heat_flux" SoilCanopyModel p.canopy.energy.shf
+@land_compute "canopy_aerodynamic_resistance" SoilCanopyModel p.canopy.energy.r_ae
+@land_compute "canopy_temperature" SoilCanopyModel Y.canopy.energy.T
+@land_compute "soil_ice" SoilCanopyModel top_center_to_surface(Y.soil.θ_i)
+
+
+#=
+
 # stored in p
 
 function compute_soil_net_radiation!(out, Y, p, t, land_model::SoilCanopyModel)
@@ -86,13 +164,7 @@ function compute_soil_water_liquid!(
     end
 end
 
-function compute_infiltration(
-    out,
-    Y,
-    p,
-    t,
-    land_model::Union{SoilCanopyModel, EnergyHydrology},
-)
+function compute_infiltration!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.soil.infiltration)
     else
@@ -100,7 +172,7 @@ function compute_infiltration(
     end
 end
 
-function compute_soilco2_diffusivity(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_soilco2_diffusivity!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.soilco2.D) # NOTE: we will need a method to compute surface co2 efflux
     else
@@ -108,7 +180,7 @@ function compute_soilco2_diffusivity(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_soilco2_source_microbe(
+function compute_soilco2_source_microbe!(
     out,
     Y,
     p,
@@ -122,7 +194,7 @@ function compute_soilco2_source_microbe(
     end
 end
 
-function compute_stomatal_conductance(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_stomatal_conductance!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.conductance.gs) # doublecheck: stomata, not canopy
     else
@@ -130,7 +202,7 @@ function compute_stomatal_conductance(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_medlyn_term(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_medlyn_term!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.conductance.medlyn_term)
     else
@@ -138,7 +210,7 @@ function compute_medlyn_term(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_canopy_transpiration(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_canopy_transpiration!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.transpiration) # doublecheck: canopy, not leaf
     else
@@ -146,7 +218,7 @@ function compute_canopy_transpiration(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_rainfall(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_rainfall!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.P_liq) # I guess this is read and put it p. not computed. curious if we should handle this differently.
     else
@@ -154,7 +226,7 @@ function compute_rainfall(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_snowfall(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_snowfall!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.P_snow) # following comment above, we could have a default getting only model output, and one also getting some inputs like drivers
     else
@@ -162,7 +234,7 @@ function compute_snowfall(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_pressure(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_pressure!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.P) # not sure if precip or pressure
     else
@@ -170,7 +242,7 @@ function compute_pressure(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_wind_speed(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_wind_speed!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.u)
     else
@@ -178,7 +250,7 @@ function compute_wind_speed(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_specific_humidity(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_specific_humidity!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.q) # check if this is correct. Also, check if Bucket has it and same name or not.
     else
@@ -186,7 +258,7 @@ function compute_specific_humidity(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_air_co2(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_air_co2!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.c_co2)
     else
@@ -194,7 +266,7 @@ function compute_air_co2(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_radiation_shortwave_down(
+function compute_radiation_shortwave_down!(
     out,
     Y,
     p,
@@ -208,7 +280,7 @@ function compute_radiation_shortwave_down(
     end
 end
 
-function compute_radiation_longwave_down(
+function compute_radiation_longwave_down!(
     out,
     Y,
     p,
@@ -222,7 +294,7 @@ function compute_radiation_longwave_down(
     end
 end
 
-function compute_photosynthesis_net_leaf(
+function compute_photosynthesis_net_leaf!(
     out,
     Y,
     p,
@@ -250,7 +322,7 @@ function compute_photosynthesis_net_canopy!(
     end
 end
 
-function compute_respiration_leaf(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_respiration_leaf!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.photosynthesis.Rd)
     else
@@ -258,7 +330,7 @@ function compute_respiration_leaf(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_vcmax25(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_vcmax25!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.photosynthesis.vcmax25)
     else
@@ -266,7 +338,7 @@ function compute_vcmax25(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_photosynthetically_active_radiation(
+function compute_photosynthetically_active_radiation!(
     out,
     Y,
     p,
@@ -280,7 +352,7 @@ function compute_photosynthetically_active_radiation(
     end
 end
 
-function compute_photosynthetically_active_radiation_absorbed(
+function compute_photosynthetically_active_radiation_absorbed!(
     out,
     Y,
     p,
@@ -294,7 +366,7 @@ function compute_photosynthetically_active_radiation_absorbed(
     end
 end
 
-function compute_photosynthetically_active_radiation_reflected(
+function compute_photosynthetically_active_radiation_reflected!(
     out,
     Y,
     p,
@@ -308,7 +380,7 @@ function compute_photosynthetically_active_radiation_reflected(
     end
 end
 
-function compute_photosynthetically_active_radiation_transmitted(
+function compute_photosynthetically_active_radiation_transmitted!(
     out,
     Y,
     p,
@@ -322,7 +394,7 @@ function compute_photosynthetically_active_radiation_transmitted(
     end
 end
 
-function compute_near_infrared_radiation(
+function compute_near_infrared_radiation!(
     out,
     Y,
     p,
@@ -336,7 +408,7 @@ function compute_near_infrared_radiation(
     end
 end
 
-function compute_near_infrared_radiation_absorbed(
+function compute_near_infrared_radiation_absorbed!(
     out,
     Y,
     p,
@@ -350,7 +422,7 @@ function compute_near_infrared_radiation_absorbed(
     end
 end
 
-function compute_near_infrared_radiation_reflected(
+function compute_near_infrared_radiation_reflected!(
     out,
     Y,
     p,
@@ -364,7 +436,7 @@ function compute_near_infrared_radiation_reflected(
     end
 end
 
-function compute_near_infrared_radiation_transmitted(
+function compute_near_infrared_radiation_transmitted!(
     out,
     Y,
     p,
@@ -378,7 +450,7 @@ function compute_near_infrared_radiation_transmitted(
     end
 end
 
-function compute_radiation_shortwave_net(
+function compute_radiation_shortwave_net!(
     out,
     Y,
     p,
@@ -392,7 +464,7 @@ function compute_radiation_shortwave_net(
     end
 end
 
-function compute_radiation_longwave_net(
+function compute_radiation_longwave_net!(
     out,
     Y,
     p,
@@ -406,7 +478,7 @@ function compute_radiation_longwave_net(
     end
 end
 
-function compute_autotrophic_respiration(
+function compute_autotrophic_respiration!(
     out,
     Y,
     p,
@@ -420,7 +492,7 @@ function compute_autotrophic_respiration(
     end
 end
 
-function compute_soilco2(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_soilco2!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(Y.soilco2.C)
     else
@@ -428,7 +500,7 @@ function compute_soilco2(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_heterotrophic_respiration(
+function compute_heterotrophic_respiration!(
     out,
     Y,
     p,
@@ -438,11 +510,11 @@ function compute_heterotrophic_respiration(
     if isnothing(out)
         return copy(p.soilco2.top_bc)
     else
-        p.soilco2.top_bc
+        out .= p.soilco2.top_bc
     end
 end
 
-function compute_soil_hydraulic_conductivity(
+function compute_soil_hydraulic_conductivity!(
     out,
     Y,
     p,
@@ -452,19 +524,19 @@ function compute_soil_hydraulic_conductivity(
     if isnothing(out)
         return copy(p.soil.K)
     else
-        p.soil.K
+        out .= p.soil.K
     end
 end
 
-function compute_soil_water_potential(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_soil_water_potential!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.soil.ψ)
     else
-        p.soil.ψ
+        out .= p.soil.ψ
     end
 end
 
-function compute_soil_thermal_conductivity(
+function compute_soil_thermal_conductivity!(
     out,
     Y,
     p,
@@ -474,19 +546,19 @@ function compute_soil_thermal_conductivity(
     if isnothing(out)
         return copy(p.soil.κ)
     else
-        p.soil.κ
+        out .= p.soil.κ
     end
 end
 
-function compute_solar_zenith_angle(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_solar_zenith_angle!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.θs)
     else
-        p.drivers.θs
+        out .= p.drivers.θs
     end
 end
 
-function compute_moisture_stress_factor(
+function compute_moisture_stress_factor!(
     out,
     Y,
     p,
@@ -496,11 +568,11 @@ function compute_moisture_stress_factor(
     if isnothing(out)
         return copy(p.canopy.hydraulics.β)
     else
-        p.canopy.hydraulics.β
+        out .= p.canopy.hydraulics.β
     end
 end
 
-function compute_canopy_water_potential(
+function compute_canopy_water_potential!(
     out,
     Y,
     p,
@@ -510,23 +582,23 @@ function compute_canopy_water_potential(
     if isnothing(out)
         return copy(p.canopy.hydraulics.ψ)
     else
-        p.canopy.hydraulics.ψ
+        out .= p.canopy.hydraulics.ψ
     end
 end
 
-function compute_cross_section(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_cross_section!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.hydraulics.fa)
     else
-        p.canopy.hydraulics.fa
+        out .= p.canopy.hydraulics.fa
     end
 end
 
-function compute_cross_section_roots(out, Y, p, t, land_model::SoilCanopyModel)
+function compute_cross_section_roots!(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.hydraulics.fa_roots)
     else
-        p.canopy.hydraulics.fa_roots
+        out .= p.canopy.hydraulics.fa_roots
     end
 end
 
@@ -538,7 +610,7 @@ function compute_area_index!(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_canopy_latent_heat_flux(
+function compute_canopy_latent_heat_flux!(
     out,
     Y,
     p,
@@ -548,11 +620,11 @@ function compute_canopy_latent_heat_flux(
     if isnothing(out)
         return copy(p.canopy.energy.lhf)
     else
-        p.canopy.canopy.lhf
+        out .= p.canopy.canopy.lhf
     end
 end
 
-function compute_canopy_sensible_heat_flux(
+function compute_canopy_sensible_heat_flux!(
     out,
     Y,
     p,
@@ -562,11 +634,11 @@ function compute_canopy_sensible_heat_flux(
     if isnothing(out)
         return copy(p.canopy.energy.shf)
     else
-        p.canopy.canopy.shf
+        out .= p.canopy.canopy.shf
     end
 end
 
-function compute_canopy_aerodynamic_resistance(
+function compute_canopy_aerodynamic_resistance!(
     out,
     Y,
     p,
@@ -576,7 +648,7 @@ function compute_canopy_aerodynamic_resistance(
     if isnothing(out)
         return copy(p.canopy.energy.r_ae)
     else
-        p.canopy.canopy.r_ae
+        out .= p.canopy.canopy.r_ae
     end
 end
 
@@ -601,3 +673,5 @@ function compute_soil_ice!(
         out .= top_center_to_surface(Y.soil.θ_i)
     end
 end
+
+=#