Merge #303

303: Varying Diffuse Radiation Fraction r=Espeer5 a=Espeer5



Co-authored-by: Edward Speer <espeer@caltech.edu>

docs/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.1"
 manifest_format = "2.0"
-project_hash = "cdeba84bbd82f27d7197cfd1d66867d60cc95b15"
+project_hash = "db0cff940cb734405cc8ed1edd362cc4e50ddc3e"
 
 [[deps.AMD]]
 deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"]

docs/Project.toml

@@ -4,6 +4,7 @@ CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
 ClimaLSM = "7884a58f-fab6-4fd0-82bb-ecfedb2d8430"
 ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"

docs/tutorials/Bucket/bucket_tutorial.jl

@@ -151,6 +151,9 @@ using ClimaLSM:
 # We also want to plot the solution
 using Plots
 
+# And we need to use the DateTime type to store reference times
+using Dates
+
 FT = Float64;
 
 # As mentioned we use CLIMAParameters for earth parameters that are
@@ -212,6 +215,11 @@ bucket_domain =
     LSMSingleColumnDomain(; zlim = (-soil_depth, 0.0), nelements = 10);
 
 
+# The PrescribedAtmosphere and PrescribedRadiation need to take in a reference 
+# time, the date of the start of the simulation. In this tutorial we will 
+# consider this January 1, 2005.
+ref_time = DateTime(2005);
+
 # To drive the system in standalone mode,
 # the user must provide
 # prescribed functions of time for the water volume flux in precipitation,
@@ -227,7 +235,8 @@ LW_d = (t) -> eltype(t)(300);
 bucket_rad = PrescribedRadiativeFluxes(
     FT,
     SW_d,
-    LW_d;
+    LW_d,
+    ref_time;
     orbital_data = Insolation.OrbitalData(),
 );
 
@@ -250,6 +259,7 @@ bucket_atmos = PrescribedAtmosphere(
     u_atmos,
     q_atmos,
     P_atmos,
+    ref_time,
     h_atmos,
 );
 

docs/tutorials/Canopy/canopy_tutorial.jl

@@ -134,7 +134,6 @@ rt_params = TwoStreamParameters{FT}(;
     λ_γ_PAR = FT(5e-7),
     λ_γ_NIR = FT(1.65e-6),
     n_layers = UInt64(20),
-    diff_perc = FT(0.5),
 )
 
 rt_model = TwoStreamModel{FT}(rt_params);

docs/tutorials/Canopy/soil_canopy_tutorial.jl

@@ -197,7 +197,6 @@ radiative_transfer_args = (;
         λ_γ_PAR = FT(5e-7),
         λ_γ_NIR = FT(1.65e-6),
         n_layers = UInt64(20),
-        diff_perc = FT(0.5),
     )
 )
 

experiments/integrated/ozark/conservation/ozark_conservation.jl

@@ -86,7 +86,6 @@ radiative_transfer_args = (;
         α_NIR_leaf = α_NIR_leaf,
         τ_NIR_leaf = τ_NIR_leaf,
         n_layers = n_layers,
-        diff_perc = diff_perc,
     )
 )
 # Set up conductance

experiments/integrated/ozark/ozark.jl

@@ -83,7 +83,6 @@ radiative_transfer_args = (;
         α_NIR_leaf = α_NIR_leaf,
         τ_NIR_leaf = τ_NIR_leaf,
         n_layers = n_layers,
-        diff_perc = diff_perc,
     )
 )
 # Set up conductance

experiments/integrated/ozark/ozark_met_drivers_FLUXNET.jl

@@ -92,14 +92,16 @@ precipitation_function(t::FT) where {FT} = p_spline(t) < 0.0 ? p_spline(t) : 0.0
 snow_precip(t) = eltype(t)(0) # this is likely not correct
 
 
-# Construct the drivers
+# Construct the drivers. For the reference time we will use the UTC time at the
+# start of the simulation
 atmos = ClimaLSM.PrescribedAtmosphere(
     precipitation_function,
     snow_precip,
     atmos_T,
     atmos_u,
     atmos_q,
     atmos_p,
+    UTC_DATETIME[1],
     atmos_h;
     c_co2 = atmos_co2,
 )
@@ -109,13 +111,14 @@ long = FT(-92.2000) # degree
 
 function zenith_angle(
     t::FT,
-    orbital_data;
+    orbital_data,
+    ref_time;
     latitude = lat,
     longitude = long,
     insol_params = earth_param_set.insol_params,
 ) where {FT}
     # This should be time in UTC
-    dt = DateTime("2005-01-01-06", "yyyy-mm-dd-HH") + Dates.Second(round(t))
+    dt = ref_time + Dates.Second(round(t))
     FT(
         instantaneous_zenith_angle(
             dt,
@@ -131,7 +134,8 @@ end
 radiation = ClimaLSM.PrescribedRadiativeFluxes(
     FT,
     SW_IN_spline,
-    LW_IN_spline;
+    LW_IN_spline,
+    UTC_DATETIME[1];
     θs = zenith_angle,
     orbital_data = Insolation.OrbitalData(),
 )

experiments/standalone/Soil/evaporation.jl

@@ -8,6 +8,7 @@ import CLIMAParameters as CP
 using RootSolvers
 using SurfaceFluxes
 using StaticArrays
+using Dates
 
 using ClimaLSM
 using ClimaLSM.Domains: Column
@@ -54,12 +55,14 @@ NIR_albedo = FT(0.4)
 z_0m = 1e-4
 z_0b = 1e-5
 
+ref_time = DateTime(2005)
 SW_d = (t) -> eltype(t)(0)
 LW_d = (t) -> eltype(t)(301.15^4 * 5.67e-8)
 radiation = PrescribedRadiativeFluxes(
     FT,
     SW_d,
-    LW_d;
+    LW_d,
+    ref_time;
     orbital_data = Insolation.OrbitalData(),
 )
 # Atmos
@@ -86,6 +89,7 @@ atmos = PrescribedAtmosphere(
     u_atmos,
     q_atmos,
     P_atmos,
+    ref_time,
     h_atmos;
     gustiness = gustiness,
 )

src/shared_utilities/drivers.jl

@@ -34,7 +34,7 @@ An abstract type of radiative drivers of land models.
 abstract type AbstractRadiativeDrivers{FT <: AbstractFloat} end
 
 """
-    PrescribedAtmosphere{FT, LP, SP, TA, UA, QA, RA, CA} <: AbstractAtmosphericDrivers{FT}
+    PrescribedAtmosphere{FT, LP, SP, TA, UA, QA, RA, CA, DT} <: AbstractAtmosphericDrivers{FT}
 
 Container for holding prescribed atmospheric drivers and other
 information needed for computing turbulent surface fluxes when
@@ -44,7 +44,7 @@ Since not all models require co2 concentration, the default for that
 is `nothing`.
 $(DocStringExtensions.FIELDS)
 """
-struct PrescribedAtmosphere{FT, LP, SP, TA, UA, QA, RA, CA} <:
+struct PrescribedAtmosphere{FT, LP, SP, TA, UA, QA, RA, CA, DT} <:
        AbstractAtmosphericDrivers{FT}
     "Precipitation (m/s) function of time: positive by definition"
     liquid_precip::LP
@@ -60,6 +60,8 @@ struct PrescribedAtmosphere{FT, LP, SP, TA, UA, QA, RA, CA} <:
     P::RA
     "CO2 concentration in atmosphere (mol/mol)"
     c_co2::CA
+    "Reference time - the datetime corresponding to t=0 for the simulation"
+    ref_time::DT
     "Reference height (m), relative to surface elevation"
     h::FT
     "Minimum wind speed (gustiness; m/s)"
@@ -71,11 +73,12 @@ struct PrescribedAtmosphere{FT, LP, SP, TA, UA, QA, RA, CA} <:
         u,
         q,
         P,
+        ref_time,
         h::FT;
         gustiness = FT(1),
         c_co2 = nothing,
     ) where {FT}
-        args = (liquid_precip, snow_precip, T, u, q, P, c_co2)
+        args = (liquid_precip, snow_precip, T, u, q, P, c_co2, ref_time)
         return new{typeof(h), typeof.(args)...}(args..., h, gustiness)
     end
 
@@ -241,29 +244,32 @@ function surface_fluxes_at_a_point(
 end
 
 """
-    PrescribedRadiativeFluxes{FT, SW, LW, T} <: AbstractRadiativeDrivers{FT}
+    PrescribedRadiativeFluxes{FT, SW, LW, DT, T, OD} <: AbstractRadiativeDrivers{FT}
 
 Container for the prescribed radiation functions needed to drive land models in standalone mode.
 $(DocStringExtensions.FIELDS)
 """
-struct PrescribedRadiativeFluxes{FT, SW, LW, T, OD} <:
+struct PrescribedRadiativeFluxes{FT, SW, LW, DT, T, OD} <:
        AbstractRadiativeDrivers{FT}
     "Downward shortwave radiation function of time (W/m^2): positive indicates towards surface"
     SW_d::SW
     "Downward longwave radiation function of time (W/m^2): positive indicates towards surface"
     LW_d::LW
+    "Reference time - the datetime corresponding to t=0 for the simulation"
+    ref_time::DT
     "Sun zenith angle, in radians"
     θs::T
     "Orbital Data for Insolation.jl"
     orbital_data::OD
     function PrescribedRadiativeFluxes(
         FT,
         SW_d,
-        LW_d;
+        LW_d,
+        ref_time;
         θs = nothing,
         orbital_data,
     )
-        args = (SW_d, LW_d, θs, orbital_data)
+        args = (SW_d, LW_d, ref_time, θs, orbital_data)
         @assert !isnothing(orbital_data)
         return new{FT, typeof.(args)...}(args...)
     end

src/standalone/Vegetation/Canopy.jl

@@ -36,6 +36,7 @@ include("./stomatalconductance.jl")
 include("./photosynthesis.jl")
 include("./radiation.jl")
 include("./canopy_parameterizations.jl")
+using Dates
 
 """
     SharedCanopyParameters{FT <: AbstractFloat, PSE}
@@ -414,7 +415,8 @@ function ClimaLSM.make_update_aux(
         (; sc, pc) = canopy.photosynthesis.parameters
 
         # Current atmospheric conditions
-        θs::FT = canopy.radiation.θs(t, canopy.radiation.orbital_data)
+        ref_time = canopy.atmos.ref_time
+        θs::FT = canopy.radiation.θs(t, canopy.radiation.orbital_data, ref_time)
         c_co2::FT = canopy.atmos.c_co2(t)
         P::FT = canopy.atmos.P(t)
         u::FT = canopy.atmos.u(t)
@@ -427,6 +429,21 @@ function ClimaLSM.make_update_aux(
         K = extinction_coeff(ld, θs)
         PAR .= compute_PAR(RT, canopy.radiation, t)
         NIR .= compute_NIR(RT, canopy.radiation, t)
+        e_sat =
+            Thermodynamics.saturation_vapor_pressure.(
+                Ref(thermo_params),
+                T,
+                Ref(Thermodynamics.Liquid()),
+            )
+        e =
+            Thermodynamics.partial_pressure_vapor.(
+                Ref(thermo_params),
+                P,
+                Ref(PhasePartition(q)),
+            )
+        rel_hum = e / e_sat
+        DOY = Dates.dayofyear(ref_time + Dates.Second(floor(Int64, t)))
+        frac_diff = @. diffuse_fraction(DOY, T, PAR + NIR, rel_hum, θs)
         APAR, ANIR = compute_absorbances(
             RT,
             PAR ./ (energy_per_photon_PAR * N_a),
@@ -436,6 +453,7 @@ function ClimaLSM.make_update_aux(
             θs,
             ground_albedo_PAR(canopy),
             ground_albedo_NIR(canopy),
+            frac_diff,
         )
 
         # update plant hydraulics aux