broadcastable solar_radiation

src/landscape.jl

@@ -1,11 +1,12 @@
-@kwdef struct SoilEnergyInputs{F,B,SP,D<:Vector{<:Number},H<:Vector{<:Number},ST,MT,EI,SW}
+@kwdef struct SoilEnergyInputs{F,B,SP,D<:Vector{<:Number},H<:Vector{<:Number},A,MT,EI,SW}
     forcing::F
     buffers::B
     soil_thermal_model::SP
     depths::D
     heights::H
+    slope::Quantity
+    albedo::A
     nodes::Vector{Float64}
-    solar_terrain::ST #TODO make just one terrain
     micro_terrain::MT
     environment_instant::EI
     runmoist::Bool
@@ -14,7 +15,7 @@ end
 
 @kwdef struct MicroForcing{
     S<:AbstractInterpolation,ZE<:AbstractInterpolation,ZS<:AbstractInterpolation,T<:AbstractInterpolation,
-    V<:AbstractInterpolation,RH<:AbstractInterpolation,CL<:AbstractInterpolation,P<:AbstractInterpolation,
+    V<:AbstractInterpolation,RH<:AbstractInterpolation,CL<:AbstractInterpolation,P<:AbstractInterpolation
 }
     interpolate_solar::S
     interpolate_zenith::ZE

src/radiation.jl

@@ -90,8 +90,8 @@ Maxwell, E. L., "A Quasi-Physical Model for Converting Hourly
 function cloud_adjust_radiation(output, cloud::AbstractArray, D_cs, B_cs, zenith::AbstractArray, doy; 
     a=0.36, b=0.64, gamma=1.0,
 )
-    (; global_total, diffuse_total, direct_total) = output.solar_radiation
-    G, D, B = (global_total, diffuse_total, direct_total)
+    (; global_horizontal, diffuse_horizontal, direct_horizontal) = output.solar_radiation
+    G, D, B = (global_horizontal, diffuse_horizontal, direct_horizontal)
     # Solar geometry
     cosz     = cos.(zenith)
     cosz_pos = max.(cosz, 0.0)

src/simulation.jl

@@ -74,8 +74,11 @@ Returns a named tuple containing:
     # TODO: this should be mandatory with no default??
     latitude
     # Objects defined above
-    solar_model = SolarProblem()
+    slope
+    aspect
+    albedo
     solar_terrain
+    solar_model = SolarProblem()
     micro_terrain
     soil_moisture_model
     soil_thermal_model
@@ -103,7 +106,7 @@ function example_microclimate_problem(;
     kw...
 )
     MicroProblem(; 
-         latitude, micro_terrain, solar_terrain, soil_moisture_model, soil_thermal_model, environment_minmax, environment_daily, kw...
+         latitude, micro_terrain, slope, aspect, albedo, solar_terrain, soil_moisture_model, soil_thermal_model, environment_minmax, environment_daily, kw...
     )
 end
 
@@ -229,9 +232,11 @@ function solve(mp::MicroProblem)
 end
 
 function solve_solar(mp::MicroProblem)
-    (; solar_model, days, hours, latitude, solar_terrain) = mp
+    (; micro_terrain, solar_model, days, hours, latitude, slope, aspect, albedo, environment_hourly, solar_terrain) = mp
+    P_atmos = environment_hourly.pressure[1] # TODO constant for now, ideally broadcast it across pressure
+    elevation = micro_terrain.elevation
     # compute clear sky solar radiation
-    solar_radiation_out = solar_radiation(solar_model; days, hours, latitude, solar_terrain)
+    solar_radiation_out = solar_radiation(solar_model, latitude, elevation, slope, aspect, albedo, P_atmos; days, hours, solar_terrain=solar_terrain)
     # limit max zenith angles to 90°
     solar_radiation_out.zenith_angle[solar_radiation_out.zenith_angle .> 90u"°"] .= 90u"°"
     solar_radiation_out.zenith_slope_angle[solar_radiation_out.zenith_slope_angle .> 90u"°"] .= 90u"°"
@@ -268,23 +273,24 @@ function adjust_for_cloud_cover(output, solar_radiation_out, days, hours)
     # adjust for cloud using Angstrom formula (formula 5.33 on P. 177 of "Climate Data and Resources" by Edward Linacre 1992
     day_of_year = repeat(days, inner=length(hours))
     zenith_angle = solar_radiation_out.zenith_angle
-    direct_total = solar_radiation_out.direct_total
-    diffuse_total = solar_radiation_out.diffuse_total
+    direct_horizontal = solar_radiation_out.direct_horizontal
+    diffuse_horizontal = solar_radiation_out.diffuse_horizontal
     cloud = output.cloud_cover
-    return (; global_solar, diffuse_fraction) = cloud_adjust_radiation(output, cloud, diffuse_total, direct_total, zenith_angle, day_of_year)
+    return (; global_solar, diffuse_fraction) = cloud_adjust_radiation(output, cloud, diffuse_horizontal, direct_horizontal, zenith_angle, day_of_year)
 end
 
 # Solves soil temperature and moisture
 function solve_soil!(output::MicroResult, mp::MicroProblem, solar_radiation_out; 
     days, hours, depths, heights
 )
-    (; solar_terrain, micro_terrain, soil_thermal_model, soil_moisture_model, environment_minmax, environment_daily, daily, initial_soil_temperature, initial_soil_moisture, runmoist, hourly_rainfall) = mp
+    (; micro_terrain, soil_thermal_model, soil_moisture_model, environment_minmax, environment_daily, environment_hourly, daily, initial_soil_temperature, initial_soil_moisture, runmoist, hourly_rainfall, slope, albedo) = mp
     (; moist_step, Campbells_b_parameter, soil_bulk_density2, soil_mineral_density2, air_entry_water_potential) = soil_moisture_model
-
     ndays = length(days)
     nhours = length(hours)
     numnodes_a = length(depths) # number of soil nodes for temperature calcs and final output
     numnodes_b = numnodes_a * 2 - 2 # number of soil nodes for soil moisture calcs
+    P_atmos = environment_hourly.pressure[1] # TODO constant for now, ideally broadcast it across pressure
+    elevation = micro_terrain.elevation
 
     # initial conditions
     if !daily && isnothing(initial_soil_temperature)
@@ -373,7 +379,7 @@ function solve_soil!(output::MicroResult, mp::MicroProblem, solar_radiation_out;
                 if i == 1 # make first hour of day equal last hour of previous iteration
                     if niter > 1
                         ∑phase .= 0.0u"J" # TODO decide whether this should happen and fix in Fortran
-                        inputs = SoilEnergyInputs(; forcing, buffers, soil_thermal_model, depths, heights, solar_terrain, micro_terrain, runmoist, nodes, environment_instant, soil_wetness)
+                        inputs = SoilEnergyInputs(; forcing, buffers, soil_thermal_model, depths, heights, slope, albedo, micro_terrain, runmoist, nodes, environment_instant, soil_wetness)
                         soiltemps = get_soil_temp_timeline(T0, inputs, i + 1)
                         T0 = soiltemps[2]
                         if iter == niter # TODO this should happen every time but at present it doesn't in Fortran version
@@ -383,7 +389,7 @@ function solve_soil!(output::MicroResult, mp::MicroProblem, solar_radiation_out;
                         end
                     end
                 else
-                    inputs = SoilEnergyInputs(; forcing, buffers, soil_thermal_model, depths, heights, solar_terrain, micro_terrain, runmoist, nodes, environment_instant, soil_wetness)
+                    inputs = SoilEnergyInputs(; forcing, buffers, soil_thermal_model, depths, heights, slope, albedo, micro_terrain, runmoist, nodes, environment_instant, soil_wetness)
                     soiltemps = get_soil_temp_timeline(T0, inputs, i)
                     T0 = soiltemps[2]
                     if iter == niter # TODO this should happen every time but at present it doesn't in Fortran version

src/soil_balance.jl

@@ -15,13 +15,12 @@ function soil_energy_balance(
     #T_K = T .* u"K"  # convert Float64 time back to unitful
     #dT_K = dT .* 60 .* u"K/minute"  # convert Float64 time back to unitful
     # extract prameters
-    (; soil_thermal_model, forcing, buffers, heights, depths, nodes, environment_instant, solar_terrain, micro_terrain, soil_wetness, runmoist) = p
+    (; soil_thermal_model, slope, albedo, forcing, buffers, heights, depths, nodes, environment_instant, micro_terrain, soil_wetness, runmoist) = p
     (; depp, wc, c) = buffers.soil_energy_balance
     (; soil_moisture, shade) = environment_instant
     # Get environmental data at time t
     (; P_atmos, tair, vel, zenr, solr, cloud, rh, zslr) = interpolate_forcings(forcing, t)
     (; roughness_height, karman_constant, dyer_constant) = micro_terrain
-    (; albedo, slope) = solar_terrain
 
     reference_height = last(heights)
     sabnew = 1.0 - albedo

test/micro_testrun_daily.jl

@@ -35,6 +35,11 @@ depths = ((DataFrame(CSV.File("$testdir/data/init_daily/DEP.csv"))[:, 2]) / 100.
 heights = [microinput[:Usrhyt], microinput[:Refhyt]]u"m" # air nodes for temperature, wind speed and humidity profile
 days2do = 30
 hours2do = days2do * 24
+slope = (microinput[:slope])*1.0u"°"
+aspect = (microinput[:azmuth])*1.0u"°"
+elevation = (microinput[:ALTT])*1.0u"m"
+albedo = (DataFrame(CSV.File("$testdir/data/init_daily/REFLS.csv"))[1, 2] * 1.0)
+P_atmos = atmospheric_pressure((microinput[:ALTT])*1.0u"m")
 
 #TODO make one terrain object via BiophysicalEcologyBase or Habitat
 micro_terrain = MicroTerrain(;
@@ -46,12 +51,7 @@ micro_terrain = MicroTerrain(;
 )
 
 solar_terrain = SolarTerrain(;
-    slope = (microinput[:slope])*1.0u"°",
-    aspect = (microinput[:azmuth])*1.0u"°",
-    elevation = (microinput[:ALTT])*1.0u"m",
     horizon_angles = (DataFrame(CSV.File("$testdir/data/init_daily/hori.csv"))[:, 2])*1.0u"°",
-    albedo = (DataFrame(CSV.File("$testdir/data/init_daily/REFLS.csv"))[1, 2] * 1.0),
-    P_atmos = atmospheric_pressure((microinput[:ALTT])*1.0u"m"),
 )
 
 soil_thermal_model = CampbelldeVriesSoilThermal(;
@@ -134,8 +134,11 @@ problem = MicroProblem(;
     depths, # soil nodes - keep spacing close near the surface
     heights, # air nodes for temperature, wind speed and humidity profile
     # Objects defined above
-    solar_model,
+    slope,
+    aspect,
+    albedo,
     solar_terrain,
+    solar_model,
     micro_terrain, #TODO combine terrains via a generic terrain in BiophysicalEcologyBase
     soil_moisture_model,
     soil_thermal_model,

test/micro_testrun_monthly.jl

@@ -33,9 +33,13 @@ LAIs = fill(0.1, length(days))
 depths = ((DataFrame(CSV.File("$testdir/data/init_monthly/DEP.csv"))[:, 2]) / 100.0)u"m"
 heights = [microinput[:Usrhyt], microinput[:Refhyt]]u"m" # air nodes for temperature, wind speed and humidity profile
 days2do = 1:12
+slope = (microinput[:slope])*1.0u"°"
+aspect = (microinput[:azmuth])*1.0u"°"
+elevation = (microinput[:ALTT])*1.0u"m"
+albedo = (DataFrame(CSV.File("$testdir/data/init_monthly/REFLS.csv"))[1, 2] * 1.0)
+P_atmos = atmospheric_pressure((microinput[:ALTT])*1.0u"m")
 
 #TODO make one terrain object via BiophysicalEcologyBase or Habitat
-#TODO make P_atmos time a varying input
 micro_terrain = MicroTerrain(;
     elevation = microinput[:ALTT] * 1.0u"m", # elevation (m)
     roughness_height = microinput[:RUF] * 1.0u"m", # roughness height for standard mode TODO dispatch based on roughness pars
@@ -45,12 +49,7 @@ micro_terrain = MicroTerrain(;
 )
 
 solar_terrain = SolarTerrain(;
-    slope = (microinput[:slope])*1.0u"°",
-    aspect = (microinput[:azmuth])*1.0u"°",
-    elevation = (microinput[:ALTT])*1.0u"m",
     horizon_angles = (DataFrame(CSV.File("$testdir/data/init_monthly/hori.csv"))[:, 2])*1.0u"°",
-    albedo = (DataFrame(CSV.File("$testdir/data/init_monthly/REFLS.csv"))[1, 2] * 1.0),
-    P_atmos = atmospheric_pressure((microinput[:ALTT])*1.0u"m"),
 )
 
 mineral_density = (CSV.File("$testdir/data/init_monthly/soilprop.csv")[1, 1][6]) * 1.0u"Mg/m^3" # soil minerals density (Mg/m3)
@@ -115,8 +114,11 @@ problem = MicroProblem(;
     depths,
     heights, # air nodes for temperature, wind speed and humidity profile
     # Objects defined above
-    solar_model,
+    slope,
+    aspect,
+    albedo,
     solar_terrain,
+    solar_model,
     micro_terrain, #TODO combine terrains via a generic terrain in BiophysicalEcologyBase
     soil_moisture_model,
     soil_thermal_model,