From accec840ab6ab887322cb8c000230daf70dbec73 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Tue, 20 Aug 2024 17:20:09 -0700 Subject: [PATCH] Implement new evaporation scheme --- .buildkite/longruns_gpu/pipeline.yml | 4 +- docs/src/APIs/Soil.md | 8 - docs/tutorials/standalone/Soil/evaporation.jl | 6 +- .../Soil/evaporation_gilat_loess.jl | 14 +- docs/tutorials/standalone/Soil/sublimation.jl | 12 +- experiments/benchmarks/bucket.jl | 2 +- experiments/integrated/fluxnet/ozark_pft.jl | 6 +- experiments/integrated/fluxnet/run_fluxnet.jl | 6 +- .../conservation/ozark_conservation.jl | 6 +- experiments/long_runs/land.jl | 12 +- experiments/long_runs/soil.jl | 3 +- .../utilities/climaland_output_dataframe.jl | 2 +- .../src/utilities/makie_plots.jl | 2 +- src/diagnostics/land_compute_methods.jl | 2 +- src/integrated/soil_canopy_model.jl | 5 +- src/standalone/Soil/Soil.jl | 5 +- src/standalone/Soil/boundary_conditions.jl | 15 +- src/standalone/Soil/energy_hydrology.jl | 439 +++++++++++++++--- .../Soil/soil_hydrology_parameterizations.jl | 84 ---- .../Vegetation/canopy_boundary_fluxes.jl | 9 +- test/standalone/Soil/climate_drivers.jl | 104 ++--- .../standalone/Soil/soil_parameterizations.jl | 20 - test/standalone/Vegetation/canopy_model.jl | 16 +- 23 files changed, 470 insertions(+), 312 deletions(-) diff --git a/.buildkite/longruns_gpu/pipeline.yml b/.buildkite/longruns_gpu/pipeline.yml index 8b0cd40b72..e57a0641ed 100644 --- a/.buildkite/longruns_gpu/pipeline.yml +++ b/.buildkite/longruns_gpu/pipeline.yml @@ -37,7 +37,7 @@ steps: artifact_paths: "land_longrun_gpu/*png" agents: slurm_gpus: 1 - slurm_time: 01:00:00 + slurm_time: 03:00:00 env: CLIMACOMMS_DEVICE: "CUDA" @@ -47,6 +47,6 @@ steps: artifact_paths: "soil_longrun_gpu/*png" agents: slurm_gpus: 1 - slurm_time: 01:00:00 + slurm_time: 03:00:00 env: CLIMACOMMS_DEVICE: "CUDA" diff --git a/docs/src/APIs/Soil.md b/docs/src/APIs/Soil.md index 95635f5058..5d375fd1e6 100644 --- a/docs/src/APIs/Soil.md +++ b/docs/src/APIs/Soil.md @@ -53,14 +53,6 @@ ClimaLand.Soil.phase_change_source ClimaLand.Soil.thermal_time ``` -## Soil Surface Parameterizations - -```@docs -ClimaLand.soil.soil_resistance -ClimaLand.Soil.dry_soil_layer_thickness -ClimaLand.Soil.soil_tortuosity -``` - ## Soil Runoff Types and Methods ```@docs diff --git a/docs/tutorials/standalone/Soil/evaporation.jl b/docs/tutorials/standalone/Soil/evaporation.jl index 2ead932021..ff0052d3fc 100644 --- a/docs/tutorials/standalone/Soil/evaporation.jl +++ b/docs/tutorials/standalone/Soil/evaporation.jl @@ -223,10 +223,8 @@ sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); # Extract the evaporation at each saved step evap = [ - parent( - sv.saveval[k].soil.turbulent_fluxes.vapor_flux .* - (1 .- sv.saveval[k].soil.ice_frac), - )[1] for k in 1:length(sol.t) + parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for + k in 1:length(sol.t) ] savepath = joinpath(pkgdir(ClimaLand), "docs/tutorials/standalone/Soil/") evaporation_data = diff --git a/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl b/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl index f5b666d587..08bcc86011 100644 --- a/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl +++ b/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl @@ -37,7 +37,7 @@ thermo_params = LP.thermodynamic_parameters(earth_param_set); # Parameters K_sat = FT(0.01 / 3600 / 24) -vg_n = FT(1.45) +vg_n = FT(1.55) vg_α = FT(1.5) hcm = vanGenuchten{FT}(; α = vg_α, n = vg_n) ν = FT(0.4) @@ -256,10 +256,8 @@ driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) cb = SciMLBase.CallbackSet(driver_cb, saving_cb) sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat) evap = [ - parent( - sv.saveval[k].soil.turbulent_fluxes.vapor_flux .* - (1 .- sv.saveval[k].soil.ice_frac), - )[1] for k in 1:length(sol.t) + parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for + k in 1:length(sol.t) ]; ## Repeat with no drainage (Ksat = 0, different BC), and with evaporation, in shorter domain @@ -323,10 +321,8 @@ cb = SciMLBase.CallbackSet(driver_cb, saving_cb) sol_no_drainage = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat) evap_no_drainage = [ - parent( - sv.saveval[k].soil.turbulent_fluxes.vapor_flux .* - (1 .- sv.saveval[k].soil.ice_frac), - )[1] for k in 1:length(sol.t) + parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for + k in 1:length(sol.t) ]; diff --git a/docs/tutorials/standalone/Soil/sublimation.jl b/docs/tutorials/standalone/Soil/sublimation.jl index b143883be4..f90fa3820a 100644 --- a/docs/tutorials/standalone/Soil/sublimation.jl +++ b/docs/tutorials/standalone/Soil/sublimation.jl @@ -209,16 +209,12 @@ sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); # Extract the evaporation at each saved step evap = [ - parent( - sv.saveval[k].soil.turbulent_fluxes.vapor_flux .* - (1 .- sv.saveval[k].soil.ice_frac), - )[1] for k in 1:length(sol.t) + parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for + k in 1:length(sol.t) ] sub = [ - parent( - sv.saveval[k].soil.turbulent_fluxes.vapor_flux .* - sv.saveval[k].soil.ice_frac, - )[1] for k in 1:length(sol.t) + parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice)[1] for + k in 1:length(sol.t) ] savepath = joinpath(pkgdir(ClimaLand), "docs/tutorials/standalone/Soil/") diff --git a/experiments/benchmarks/bucket.jl b/experiments/benchmarks/bucket.jl index c62fc7d96a..2bb4f933ac 100644 --- a/experiments/benchmarks/bucket.jl +++ b/experiments/benchmarks/bucket.jl @@ -224,7 +224,7 @@ if ClimaComms.device() isa ClimaComms.CUDADevice end if get(ENV, "BUILDKITE_PIPELINE_SLUG", nothing) == "climaland-benchmark" - PREVIOUS_BEST_TIME = 3.6 + PREVIOUS_BEST_TIME = 2.8 if average_timing_s > PREVIOUS_BEST_TIME + std_timing_s @info "Possible performance regression, previous average time was $(PREVIOUS_BEST_TIME)" elseif average_timing_s < PREVIOUS_BEST_TIME - std_timing_s diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl index 92f5ba968b..4a3cd83d89 100644 --- a/experiments/integrated/fluxnet/ozark_pft.jl +++ b/experiments/integrated/fluxnet/ozark_pft.jl @@ -455,8 +455,10 @@ T = ] .* (1e3 * 24 * 3600) E = [ - parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux)[1] for - k in 1:length(sol.t) + parent( + sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq .+ + sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice, + )[1] for k in 1:length(sol.t) ] .* (1e3 * 24 * 3600) ET_model = T .+ E if drivers.LE.status == absent diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl index caf0cd05f0..4ceed3ca3d 100644 --- a/experiments/integrated/fluxnet/run_fluxnet.jl +++ b/experiments/integrated/fluxnet/run_fluxnet.jl @@ -399,8 +399,10 @@ T = ] .* (1e3 * 24 * 3600) E = [ - parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux)[1] for - k in 1:length(sol.t) + parent( + sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq .+ + sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice, + )[1] for k in 1:length(sol.t) ] .* (1e3 * 24 * 3600) ET_model = T .+ E if drivers.LE.status == absent diff --git a/experiments/integrated/performance/conservation/ozark_conservation.jl b/experiments/integrated/performance/conservation/ozark_conservation.jl index 1b4e6d7150..91fcc08cee 100644 --- a/experiments/integrated/performance/conservation/ozark_conservation.jl +++ b/experiments/integrated/performance/conservation/ozark_conservation.jl @@ -104,8 +104,10 @@ for float_type in (Float32, Float64) # Evaporation E = [ - parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux)[1] for - k in 2:length(sol.t) + parent( + sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq .+ + sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice, + )[1] for k in 2:length(sol.t) ] # Root sink term: a positive root extraction is a sink term for soil; add minus sign root_sink = diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl index 8aa50421a8..760c1dadcd 100644 --- a/experiments/long_runs/land.jl +++ b/experiments/long_runs/land.jl @@ -613,7 +613,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15)) t0, ref_time; output_writer = nc_writer, - output_vars = :long, + output_vars = :short, ) diagnostic_handler = @@ -628,7 +628,7 @@ end function setup_and_solve_problem(; greet = false) t0 = 0.0 - tf = 60 * 60.0 * 24 * 7 # keep short until it runs! * 365 + tf = 60 * 60.0 * 24 * 14 Δt = 900.0 nelements = (101, 15) if greet @@ -641,12 +641,12 @@ function setup_and_solve_problem(; greet = false) prob, cb = setup_prob(t0, tf, Δt; nelements) # Define timestepper and ODE algorithm - stepper = CTS.ARS343() + stepper = CTS.ARS111() ode_algo = CTS.IMEXAlgorithm( stepper, CTS.NewtonsMethod( - max_iters = 1, - update_j = CTS.UpdateEvery(CTS.NewTimeStep), + max_iters = 3, + update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), ), ) SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb, adaptive = false) @@ -657,7 +657,7 @@ setup_and_solve_problem(; greet = true); # read in diagnostics and make some plots! #### ClimaAnalysis #### simdir = ClimaAnalysis.SimDir(outdir) -short_names = ["gpp", "swc", "si", "sie"] +short_names = ["gpp", "ct", "lai", "swc", "si"] for short_name in short_names var = get(simdir; short_name) times = ClimaAnalysis.times(var) diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl index f0550f9499..66f22c3acd 100644 --- a/experiments/long_runs/soil.jl +++ b/experiments/long_runs/soil.jl @@ -425,6 +425,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15)) t0, ref_time; output_writer = nc_writer, + average_period = :monthly, ) diagnostic_handler = @@ -439,7 +440,7 @@ end function setup_and_solve_problem(; greet = false) t0 = 0.0 - tf = 60 * 60.0 * 24 * 60 # keep short until it runs! * 365 + tf = 60 * 60.0 * 24 * 340 Δt = 900.0 nelements = (101, 15) if greet diff --git a/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl b/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl index 26323a00a8..166c20f03a 100644 --- a/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl +++ b/lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl @@ -56,7 +56,7 @@ function make_output_df( (1, :soil, :turbulent_fluxes, :lhf), collect(map(i -> (i, :soil, :T), 1:20)), # 20 shouldn't be hard-coded, but an arg, equal to n layers collect(map(i -> (i, :soil, :θ_l), 1:20)), - (1, :soil, :turbulent_fluxes, :vapor_flux), + (1, :soil, :turbulent_fluxes, :vapor_flux_liq), (1, :canopy, :sif, :SIF), ) diff --git a/lib/ClimaLandSimulations/src/utilities/makie_plots.jl b/lib/ClimaLandSimulations/src/utilities/makie_plots.jl index d9cf9910ef..a6efe31fca 100644 --- a/lib/ClimaLandSimulations/src/utilities/makie_plots.jl +++ b/lib/ClimaLandSimulations/src/utilities/makie_plots.jl @@ -72,7 +72,7 @@ function timeseries_fluxes_fig( p_ET_m = lines!( ax_W, datetime2unix.(climaland.DateTime), - (climaland.vapor_flux .* 1e3 .* 24 .* 3600) .+ + (climaland.vapor_flux_liq .* 1e3 .* 24 .* 3600) .+ (climaland.transpiration .* 1e3 .* 24 .* 3600), color = :blue, ) # not sure about units diff --git a/src/diagnostics/land_compute_methods.jl b/src/diagnostics/land_compute_methods.jl index 041412640b..3532d949c8 100644 --- a/src/diagnostics/land_compute_methods.jl +++ b/src/diagnostics/land_compute_methods.jl @@ -136,7 +136,7 @@ end @diagnostic_compute "soil_aerodynamic_resistance" SoilCanopyModel p.soil.turbulent_fluxes.r_ae @diagnostic_compute "soil_latent_heat_flux" SoilCanopyModel p.soil.turbulent_fluxes.lhf @diagnostic_compute "soil_sensible_heat_flux" SoilCanopyModel p.soil.turbulent_fluxes.shf -@diagnostic_compute "vapor_flux" SoilCanopyModel p.soil.turbulent_fluxes.vapor_flux +@diagnostic_compute "vapor_flux" SoilCanopyModel p.soil.turbulent_fluxes.vapor_flux_liq # should add ice here # Soil - SoilCO2 function compute_heterotrophic_respiration!( diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl index ee49351665..b7c40688bc 100644 --- a/src/integrated/soil_canopy_model.jl +++ b/src/integrated/soil_canopy_model.jl @@ -418,11 +418,8 @@ function soil_boundary_fluxes!( bc = soil.boundary_conditions.top p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t) Soil.Runoff.update_runoff!(p, bc.runoff, Y, t, soil) - # Multiply the vapor flux by 1 - p.soil.ice_frac to get - # the approximated evaporation of liquid water @. p.soil.top_bc.water = - p.soil.infiltration + - p.soil.turbulent_fluxes.vapor_flux * (1 - p.soil.ice_frac) + p.soil.infiltration + p.soil.turbulent_fluxes.vapor_flux_liq @. p.soil.top_bc.heat = -p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf end diff --git a/src/standalone/Soil/Soil.jl b/src/standalone/Soil/Soil.jl index f3350a9ae1..80fc58bcd2 100644 --- a/src/standalone/Soil/Soil.jl +++ b/src/standalone/Soil/Soil.jl @@ -64,7 +64,9 @@ import ClimaCore.MatrixFields: @name, ⋅ import ..Parameters as LP import ClimaCore: Fields, Operators, Geometry, Spaces using Thermodynamics - +using SurfaceFluxes +using StaticArrays +import SurfaceFluxes.Parameters as SFP import ClimaLand.Domains: Column, HybridBox, SphericalShell import ClimaLand: AbstractImExModel, @@ -90,6 +92,7 @@ import ClimaLand: surface_emissivity, surface_height, surface_resistance, + turbulent_fluxes, get_drivers export RichardsModel, RichardsParameters, diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl index 1096b61dfb..c6aeac5af3 100644 --- a/src/standalone/Soil/boundary_conditions.jl +++ b/src/standalone/Soil/boundary_conditions.jl @@ -662,7 +662,6 @@ boundary_vars( ::ClimaLand.TopBoundary, ) = ( :turbulent_fluxes, - :ice_frac, :R_n, :top_bc, :sfc_scratch, @@ -692,7 +691,6 @@ boundary_var_domain_names( :surface, :surface, :surface, - :surface, Runoff.runoff_var_domain_names(bc.runoff)..., ) """ @@ -717,8 +715,10 @@ boundary_var_types( }, ::ClimaLand.TopBoundary, ) where {FT} = ( - NamedTuple{(:lhf, :shf, :vapor_flux, :r_ae), Tuple{FT, FT, FT, FT}}, - FT, + NamedTuple{ + (:lhf, :shf, :vapor_flux_liq, :r_ae, :vapor_flux_ice), + Tuple{FT, FT, FT, FT, FT}, + }, FT, NamedTuple{(:water, :heat), Tuple{FT, FT}}, FT, @@ -764,12 +764,9 @@ function soil_boundary_fluxes!( p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, model, Y, p, t) p.soil.R_n .= net_radiation(bc.radiation, model, Y, p, t) update_runoff!(p, bc.runoff, Y, t, model) - # We do not model the energy flux from infiltration. We multiply - # the vapor flux by the ice fraction in order to get the contribution - # from liquid water + # We do not model the energy flux from infiltration. @. p.soil.top_bc.water = - p.soil.infiltration + - p.soil.turbulent_fluxes.vapor_flux * (1 - p.soil.ice_frac) + p.soil.infiltration + p.soil.turbulent_fluxes.vapor_flux_liq @. p.soil.top_bc.heat = p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf end diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl index f562d596c8..8463284c8a 100644 --- a/src/standalone/Soil/energy_hydrology.jl +++ b/src/standalone/Soil/energy_hydrology.jl @@ -676,7 +676,7 @@ function ClimaLand.source!( z = model.domain.fields.z Δz_top = model.domain.fields.Δz_top # this returns the center-face distance, not layer thickness @. dY.soil.θ_i += - -p.soil.turbulent_fluxes.vapor_flux * p.soil.ice_frac * _ρ_l / _ρ_i * + -p.soil.turbulent_fluxes.vapor_flux_ice * _ρ_l / _ρ_i * heaviside(z + 2 * Δz_top) # only apply to top layer, recall that z is negative end @@ -712,23 +712,87 @@ function ClimaLand.surface_temperature( end """ - ClimaLand.surface_resistance( + ClimaLand.surface_emissivity( model::EnergyHydrology{FT}, Y, p, - t, ) where {FT} -Returns the surface resistance field of the +Returns the surface emissivity field of the `EnergyHydrology` soil model. """ -function ClimaLand.surface_resistance( +function ClimaLand.surface_emissivity( model::EnergyHydrology{FT}, Y, p, +) where {FT} + return model.parameters.emissivity +end + +""" + ClimaLand.surface_albedo( + model::EnergyHydrology{FT}, + Y, + p, + ) where {FT} + +Returns the surface albedo field of the +`EnergyHydrology` soil model. +""" +function ClimaLand.surface_albedo(model::EnergyHydrology{FT}, Y, p) where {FT} + return @. FT( + 0.5 * model.parameters.PAR_albedo + 0.5 * model.parameters.NIR_albedo, + ) +end + +""" + ClimaLand.surface_height( + model::EnergyHydrology{FT}, + Y, + p, + ) where {FT} + +Returns the surface height of the `EnergyHydrology` model. +""" +function ClimaLand.surface_height(model::EnergyHydrology{FT}, Y, p) where {FT} + return model.domain.fields.z_sfc +end + +function ClimaLand.get_drivers(model::EnergyHydrology) + bc = model.boundary_conditions.top + if typeof(bc) <: AtmosDrivenFluxBC{ + <:PrescribedAtmosphere, + <:AbstractRadiativeDrivers, + <:AbstractRunoffModel, + } + return (bc.atmos, bc.radiation) + else + return () + end +end + + + +function turbulent_fluxes( + atmos::PrescribedAtmosphere, + model::EnergyHydrology{FT}, + Y::ClimaCore.Fields.FieldVector, + p::NamedTuple, t, ) where {FT} - (; ν, θ_r, d_ds, earth_param_set, hydrology_cm) = model.parameters + # Obtain surface quantities needed for computation; these should not allocate + T_sfc = ClimaLand.surface_temperature(model, Y, p, t) + h_sfc = ClimaLand.surface_height(model, Y, p) + d_sfc = ClimaLand.displacement_height(model, Y, p) + u_air = p.drivers.u + h_air = atmos.h + (; K_sat, ν, θ_r, hydrology_cm, z_0m, z_0b, Ω, γ, γT_ref, earth_param_set) = + model.parameters + hydrology_cm_sfc = ClimaLand.Domains.top_center_to_surface(hydrology_cm) + K_sat_sfc = ClimaLand.Domains.top_center_to_surface(K_sat) + θ_i_sfc = ClimaLand.Domains.top_center_to_surface(Y.soil.θ_i) + ν_sfc = ClimaLand.Domains.top_center_to_surface(ν) + θ_r_sfc = ClimaLand.Domains.top_center_to_surface(θ_r) θ_l_sfc = p.soil.sfc_scratch ClimaLand.Domains.linear_interpolation_to_surface!( θ_l_sfc, @@ -736,56 +800,272 @@ function ClimaLand.surface_resistance( model.domain.fields.z, model.domain.fields.Δz_top, ) - # These are non-allocating - θ_i_sfc = ClimaLand.Domains.top_center_to_surface(Y.soil.θ_i) - hydrology_cm_sfc = ClimaLand.Domains.top_center_to_surface(hydrology_cm) - ν_sfc = ClimaLand.Domains.top_center_to_surface(ν) - θ_r_sfc = ClimaLand.Domains.top_center_to_surface(θ_r) - return ClimaLand.Soil.soil_resistance.( + return soil_turbulent_fluxes_at_a_point.( + T_sfc, θ_l_sfc, θ_i_sfc, + h_sfc, + d_sfc, hydrology_cm_sfc, ν_sfc, θ_r_sfc, - d_ds, - earth_param_set, + K_sat_sfc, + p.drivers.thermal_state, + u_air, + h_air, + atmos.gustiness, + z_0m, + z_0b, + Ω, + γ, + γT_ref, + Ref(earth_param_set), ) end """ - ClimaLand.surface_emissivity( - model::EnergyHydrology{FT}, - Y, - p, - ) where {FT} - -Returns the surface emissivity field of the -`EnergyHydrology` soil model. + soil_turbulent_fluxes_at_a_point( + T_sfc::FT, + θ_l_sfc::FT, + θ_i_sfc::FT, + h_sfc::FT, + d_sfc::FT, + hydrology_cm_sfc::C, + ν_sfc::FT, + θ_r_sfc::FT, + K_sat_sfc::FT, + thermal_state_air::Thermodynamics.PhaseEquil{FT}, + u_air::FT, + h_air::FT, + gustiness::FT, + z_0m::FT, + z_0b::FT, + Ω::FT, + γ::FT, + γT_ref,::FT + earth_param_set::EP + ) where {FT <: AbstractFloat, C, EP} + +Computes turbulent surface fluxes for soil at a point on a surface given +(1) Surface state conditions (`T_sfc`, `θ_l_sfc`, `θ_i_sfc`) +(2) Surface properties, such as the topographical height of the surface (`h_sfc`), + the displacement height (`d_sfc`), hydraulic parameters (`hydrology_cm_sfc`, + `ν_sfc, `θ_r_sfc`, `K_sat_sfc`) +(4) Atmospheric state conditions (`thermal_state_air`, `u_air`) +(5) Height corresponding to where atmospheric state is measured (`h_air`) +(6) Parameters: `gustiness`, roughness lengths `z_0m`, `z_0b`, several + required to compute the soil conductivity `Ω`, γ`, γT_ref`, and the + `earth_param_set` + +This returns an energy flux and a liquid water volume flux, stored in +a tuple with self explanatory keys. """ -function ClimaLand.surface_emissivity( - model::EnergyHydrology{FT}, - Y, - p, -) where {FT} - return model.parameters.emissivity +function soil_turbulent_fluxes_at_a_point( + T_sfc::FT, + θ_l_sfc::FT, + θ_i_sfc::FT, + h_sfc::FT, + d_sfc::FT, + hydrology_cm_sfc, + ν_sfc::FT, + θ_r_sfc::FT, + K_sat_sfc::FT, + thermal_state_air::Thermodynamics.PhaseEquil{FT}, + u_air::FT, + h_air::FT, + gustiness::FT, + z_0m::FT, + z_0b::FT, + Ω::FT, + γ::FT, + γT_ref::FT, + earth_param_set::P, +) where {FT <: AbstractFloat, P} + thermo_params = LP.thermodynamic_parameters(earth_param_set) + # Estimate surface air density + ρ_sfc::FT = ClimaLand.compute_ρ_sfc(thermo_params, thermal_state_air, T_sfc) + # Compute saturated specific humidities at surface over ice and liquid water + q_sat_ice::FT = Thermodynamics.q_vap_saturation_generic( + thermo_params, + T_sfc, + ρ_sfc, + Thermodynamics.Ice(), + ) + q_sat_liq::FT = Thermodynamics.q_vap_saturation_generic( + thermo_params, + T_sfc, + ρ_sfc, + Thermodynamics.Liquid(), + ) + + # Evaporation: + thermal_state_sfc::Thermodynamics.PhaseEquil{FT} = + Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sat_liq)# use to get potential evaporation E0, r_ae (weak dependence on q). + + # SurfaceFluxes.jl expects a relative difference between where u_air = 0 + # and the atmosphere height. Here, we assume h and h_sfc are measured + # relative to a common reference. Then d_sfc + h_sfc + z_0m is the apparent + # source of momentum, and + # Δh ≈ h_air - d_sfc - h_sfc is the relative height difference between the + # apparent source of momentum and the atmosphere height. + + # In this we have neglected z_0m and z_0b (i.e. assumed they are small + # compared to Δh). + state_sfc = SurfaceFluxes.StateValues( + FT(0), + SVector{2, FT}(0, 0), + thermal_state_sfc, + ) + state_air = SurfaceFluxes.StateValues( + h_air - d_sfc - h_sfc, + SVector{2, FT}(u_air, 0), + thermal_state_air, + ) + + # State containers + states = SurfaceFluxes.ValuesOnly( + state_air, + state_sfc, + z_0m, + z_0b, + gustiness = gustiness, + ) + surface_flux_params = LP.surface_fluxes_parameters(earth_param_set) + conditions = SurfaceFluxes.surface_conditions( + surface_flux_params, + states; + tol_neutral = SFP.cp_d(surface_flux_params) / 100000, + ) + _LH_v0::FT = LP.LH_v0(earth_param_set) + _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set) + cp_m::FT = Thermodynamics.cp_m(thermo_params, thermal_state_air) + T_air::FT = Thermodynamics.air_temperature(thermo_params, thermal_state_air) + ρ_air::FT = Thermodynamics.air_density(thermo_params, thermal_state_air) + q_air::FT = + Thermodynamics.total_specific_humidity(thermo_params, thermal_state_air) + + ΔT::FT = T_air - T_sfc + r_ae_liq::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states)) + + E0::FT = + SurfaceFluxes.evaporation(surface_flux_params, states, conditions.Ch)# potential evaporation rate, mass flux + Ẽ0::FT = E0 / _ρ_liq + Ẽ::FT = Ẽ0 + if q_air < q_sat_liq # adjust potential evaporation rate to account for soil resistance + K_sfc::FT = + impedance_factor(θ_i_sfc / (θ_l_sfc + θ_i_sfc - θ_r_sfc), Ω) * + viscosity_factor(T_sfc, γ, γT_ref) * + hydraulic_conductivity( + hydrology_cm_sfc, + K_sat_sfc, + effective_saturation(ν_sfc, θ_l_sfc, θ_r_sfc), + ) + K_c::FT = hydraulic_conductivity( + hydrology_cm_sfc, + K_sat_sfc, + hydrology_cm_sfc.S_c, + ) + x::FT = 4 * K_sfc * (1 + Ẽ0 / (4 * K_c)) + Ẽ *= x / (Ẽ0 + x) + else + Ẽ *= 0 # condensation, set to zero + end + H_liq::FT = -ρ_air * cp_m * ΔT / r_ae_liq + + # Sublimation: + thermal_state_sfc = + Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sat_ice) + β_i::FT = FT(1) + if q_air < q_sat_ice # sublimation, adjust β + β_i *= (θ_i_sfc / ν_sfc)^4 + else + β_i *= 0 # frost, set to zero + end + + state_sfc = SurfaceFluxes.StateValues( + FT(0), + SVector{2, FT}(0, 0), + thermal_state_sfc, + ) + + states = SurfaceFluxes.ValuesOnly( + state_air, + state_sfc, + z_0m, + z_0b, + beta = β_i, + gustiness = gustiness, + ) + conditions = SurfaceFluxes.surface_conditions( + surface_flux_params, + states; + tol_neutral = SFP.cp_d(surface_flux_params) / 100000, + ) + S::FT = + SurfaceFluxes.evaporation(surface_flux_params, states, conditions.Ch)# sublimation rate, mass flux + S̃::FT = S / _ρ_liq + + r_ae_ice::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states)) + H_ice::FT = -ρ_air * cp_m * ΔT / r_ae_ice + + # Heat fluxes + LH::FT = _LH_v0 * (Ẽ + S̃) * _ρ_liq + H::FT = (H_ice + H_liq) / 2 + r_ae::FT = 2 * r_ae_liq * r_ae_ice / (r_ae_liq + r_ae_ice) # implied by definition of H + return ( + lhf = LH, + shf = H, + vapor_flux_liq = Ẽ, + r_ae = r_ae, + vapor_flux_ice = S̃, + ) end + +# For Swenson/Lawrence 2014 resistance parameterization +#= """ - ClimaLand.surface_albedo( + ClimaLand.surface_resistance( model::EnergyHydrology{FT}, Y, p, + t, ) where {FT} -Returns the surface albedo field of the +Returns the surface resistance field of the `EnergyHydrology` soil model. """ -function ClimaLand.surface_albedo(model::EnergyHydrology{FT}, Y, p) where {FT} - return @. FT( - 0.5 * model.parameters.PAR_albedo + 0.5 * model.parameters.NIR_albedo, +function ClimaLand.surface_resistance( + model::EnergyHydrology{FT}, + Y, + p, + t, +) where {FT} + (; ν, θ_r, d_ds, earth_param_set, hydrology_cm) = model.parameters + θ_l_sfc = p.soil.sfc_scratch + ClimaLand.Domains.linear_interpolation_to_surface!( + θ_l_sfc, + p.soil.θ_l, + model.domain.fields.z, + model.domain.fields.Δz_top, + ) + # These are non-allocating + θ_i_sfc = ClimaLand.Domains.top_center_to_surface(Y.soil.θ_i) + hydrology_cm_sfc = ClimaLand.Domains.top_center_to_surface(hydrology_cm) + ν_sfc = ClimaLand.Domains.top_center_to_surface(ν) + θ_r_sfc = ClimaLand.Domains.top_center_to_surface(θ_r) + return ClimaLand.Soil.soil_resistance.( + θ_l_sfc, + θ_i_sfc, + hydrology_cm_sfc, + ν_sfc, + θ_r_sfc, + d_ds, + earth_param_set, ) end + """ ClimaLand.surface_specific_humidity( model::EnergyHydrology{FT}, @@ -858,29 +1138,84 @@ function ClimaLand.surface_specific_humidity( return @. (1 - ice_frac) * q_over_liq + ice_frac * q_over_ice end +""" + ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT)::FT where {FT} + +Computes and returns the ice fraction, which is the +fraction of the vapor flux that is due to sublimation, and +the fraction of the humidity in the air due to ice, as + +f = S_i/(S_i+S_l) +This same fraction is used to estimate the specific humidity, i.e. +q = q_over_ice * f + q_over_water * (1-f). """ - ClimaLand.surface_height( - model::EnergyHydrology{FT}, - Y, - p, - ) where {FT} +function ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT)::FT where {FT} + S_l = effective_saturation(ν, θ_l, θ_r) + S_i = effective_saturation(ν, θ_i, θ_r) + f = S_i / (S_i + S_l) + return f +end -Returns the surface height of the `EnergyHydrology` model. """ -function ClimaLand.surface_height(model::EnergyHydrology{FT}, Y, p) where {FT} - return model.domain.fields.z_sfc + soil_tortuosity(θ_l::FT, θ_i::FT, ν::FT) where {FT} + +Computes the tortuosity of water vapor in a porous medium, +as a function of porosity `ν` and the volumetric liquid water +and ice contents, `θ_l` and `θ_i`. + +See Equation (1) of : Shokri, N., P. Lehmann, and +D. Or (2008), Effects of hydrophobic layers on evaporation from +porous media, Geophys. Res. Lett., 35, L19407, doi:10.1029/ +2008GL035230. +""" +function soil_tortuosity(θ_l::FT, θ_i::FT, ν::FT) where {FT} + safe_θ_a = max(ν - θ_l - θ_i, eps(FT)) + return safe_θ_a^(FT(2.5)) / ν end -function ClimaLand.get_drivers(model::EnergyHydrology) - bc = model.boundary_conditions.top - if typeof(bc) <: AtmosDrivenFluxBC{ - <:PrescribedAtmosphere, - <:AbstractRadiativeDrivers, - <:AbstractRunoffModel, - } - return (bc.atmos, bc.radiation) - else - return () - end +""" + soil_resistance(θ_l::FT, + θ_i::FT, + hydrology_cm::C, + ν::FT, + θ_r::FT, + d_ds::FT, + earth_param_set::EP, + ) where {FT, EP, C} + +Computes the resistance of the top of the soil column to +water vapor diffusion, as a function of the surface +volumetric liquid water fraction `θ_l`, the augmented +liquid water fraction `ϑ_l`, the volumetric ice water +fraction `θ_i`, and other soil parameters. +""" +function soil_resistance( + θ_l::FT, + θ_i::FT, + hydrology_cm::C, + ν::FT, + θ_r::FT, + d_ds::FT, + earth_param_set::EP, +) where {FT, EP, C} + (; S_c) = hydrology_cm + _D_vapor = FT(LP.D_vapor(earth_param_set)) + S_w = effective_saturation(ν, θ_l + θ_i, θ_r) + τ_a = soil_tortuosity(θ_l, θ_i, ν) + dsl::FT = dry_soil_layer_thickness(S_w, S_c, d_ds) + r_soil = dsl / (_D_vapor * τ_a) # [s\m] + return r_soil +end + +""" + dry_soil_layer_thickness(S_w::FT, S_c::FT, d_ds::FT)::FT where {FT} + +Returns the maximum dry soil layer thickness that can develop under vapor flux; +this is used when computing the soil resistance to vapor flux according to +Swenson et al (2012)/Sakaguchi and Zeng (2009). +""" +function dry_soil_layer_thickness(S_w::FT, S_c::FT, d_ds::FT)::FT where {FT} + return S_w < S_c ? d_ds * (S_c - S_w) / S_c : FT(0) end +=# diff --git a/src/standalone/Soil/soil_hydrology_parameterizations.jl b/src/standalone/Soil/soil_hydrology_parameterizations.jl index b677c19bd5..30e0974bf7 100644 --- a/src/standalone/Soil/soil_hydrology_parameterizations.jl +++ b/src/standalone/Soil/soil_hydrology_parameterizations.jl @@ -7,9 +7,6 @@ export volumetric_liquid_fraction, impedance_factor, viscosity_factor, dψdϑ, - dry_soil_layer_thickness, - soil_resistance, - soil_tortuosity, is_saturated """ volumetric_liquid_fraction(ϑ_l::FT, ν_eff::FT, θ_r::FT) where {FT} @@ -297,87 +294,6 @@ function viscosity_factor(T::FT, γ::FT, γT_ref::FT) where {FT} return Theta end -""" - ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT)::FT where {FT} - -Computes and returns the ice fraction, which is the -fraction of the vapor flux that is due to sublimation, and -the fraction of the humidity in the air due to ice, as - -f = S_i/(S_i+S_l) - -This same fraction is used to estimate the specific humidity, i.e. -q = q_over_ice * f + q_over_water * (1-f). -""" -function ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT)::FT where {FT} - S_l = effective_saturation(ν, θ_l, θ_r) - S_i = effective_saturation(ν, θ_i, θ_r) - f = S_i / (S_i + S_l) - return f -end - -""" - soil_tortuosity(θ_l::FT, θ_i::FT, ν::FT) where {FT} - -Computes the tortuosity of water vapor in a porous medium, -as a function of porosity `ν` and the volumetric liquid water -and ice contents, `θ_l` and `θ_i`. - -See Equation (1) of : Shokri, N., P. Lehmann, and -D. Or (2008), Effects of hydrophobic layers on evaporation from -porous media, Geophys. Res. Lett., 35, L19407, doi:10.1029/ -2008GL035230. -""" -function soil_tortuosity(θ_l::FT, θ_i::FT, ν::FT) where {FT} - safe_θ_a = max(ν - θ_l - θ_i, eps(FT)) - return safe_θ_a^(FT(2.5)) / ν -end - -""" - soil_resistance(θ_l::FT, - θ_i::FT, - hydrology_cm::C, - ν::FT, - θ_r::FT, - d_ds::FT, - earth_param_set::EP, - ) where {FT, EP, C} - -Computes the resistance of the top of the soil column to -water vapor diffusion, as a function of the surface -volumetric liquid water fraction `θ_l`, the augmented -liquid water fraction `ϑ_l`, the volumetric ice water -fraction `θ_i`, and other soil parameters. -""" -function soil_resistance( - θ_l::FT, - θ_i::FT, - hydrology_cm::C, - ν::FT, - θ_r::FT, - d_ds::FT, - earth_param_set::EP, -) where {FT, EP, C} - (; S_c) = hydrology_cm - _D_vapor = FT(LP.D_vapor(earth_param_set)) - S_w = effective_saturation(ν, θ_l + θ_i, θ_r) - τ_a = soil_tortuosity(θ_l, θ_i, ν) - dsl::FT = dry_soil_layer_thickness(S_w, S_c, d_ds) - r_soil = dsl / (_D_vapor * τ_a) # [s\m] - return r_soil -end - -""" - dry_soil_layer_thickness(S_w::FT, S_c::FT, d_ds::FT)::FT where {FT} - -Returns the maximum dry soil layer thickness that can develop under vapor flux; -this is used when computing the soil resistance to vapor flux according to -Swenson et al (2012)/Sakaguchi and Zeng (2009). -""" -function dry_soil_layer_thickness(S_w::FT, S_c::FT, d_ds::FT)::FT where {FT} - return S_w < S_c ? d_ds * (S_c - S_w) / S_c : FT(0) -end - """ is_saturated(twc::FT, ν::FT) where {FT} diff --git a/src/standalone/Vegetation/canopy_boundary_fluxes.jl b/src/standalone/Vegetation/canopy_boundary_fluxes.jl index 9ee9518227..1e20b98c0a 100644 --- a/src/standalone/Vegetation/canopy_boundary_fluxes.jl +++ b/src/standalone/Vegetation/canopy_boundary_fluxes.jl @@ -9,7 +9,8 @@ import ClimaLand: surface_evaporative_scaling, surface_height, surface_resistance, - displacement_height + displacement_height, + turbulent_fluxes """ @@ -168,7 +169,7 @@ function canopy_boundary_fluxes!( i_end = canopy.hydraulics.n_stem + canopy.hydraulics.n_leaf # Compute transpiration, SHF, LHF - canopy_tf = canopy_turbulent_fluxes(atmos, canopy, Y, p, t) + canopy_tf = ClimaLand.turbulent_fluxes(atmos, canopy, Y, p, t) transpiration .= canopy_tf.vapor_flux shf .= canopy_tf.shf lhf .= canopy_tf.lhf @@ -213,7 +214,7 @@ function canopy_boundary_fluxes!( end """ - function canopy_turbulent_fluxes( + function ClimaLand.turbulent_fluxes( atmos::PrescribedAtmosphere, model::CanopyModel, Y::ClimaCore.Fields.FieldVector, @@ -229,7 +230,7 @@ because the canopy requires a different resistance for vapor and sensible heat fluxes, and the resistances depend on ustar, which we must compute using SurfaceFluxes before adjusting to account for these resistances. """ -function canopy_turbulent_fluxes( +function ClimaLand.turbulent_fluxes( atmos::PrescribedAtmosphere, model::CanopyModel, Y::ClimaCore.Fields.FieldVector, diff --git a/test/standalone/Soil/climate_drivers.jl b/test/standalone/Soil/climate_drivers.jl index 15ca894c8d..8f6bba7b78 100644 --- a/test/standalone/Soil/climate_drivers.jl +++ b/test/standalone/Soil/climate_drivers.jl @@ -38,8 +38,6 @@ for FT in (Float32, Float64) hcm = vanGenuchten{FT}(; α = vg_α, n = vg_n) θ_r = FT(0.1) S_c = hcm.S_c - @test Soil.dry_soil_layer_thickness(FT(1), S_c, FT(1)) == FT(0) - @test Soil.dry_soil_layer_thickness(FT(0), S_c, FT(1)) == FT(1) ν_ss_om = FT(0.0) ν_ss_quartz = FT(1.0) @@ -131,7 +129,7 @@ for FT in (Float32, Float64) :θs, ) @test propertynames(p.soil.turbulent_fluxes) == - (:lhf, :shf, :vapor_flux, :r_ae) + (:lhf, :shf, :vapor_flux_liq, :r_ae, :vapor_flux_ice) @test propertynames(p.soil) == ( :K, :ψ, @@ -139,7 +137,6 @@ for FT in (Float32, Float64) :T, :κ, :turbulent_fluxes, - :ice_frac, :R_n, :top_bc, :sfc_scratch, @@ -196,57 +193,11 @@ for FT in (Float32, Float64) ) T_sfc = ClimaCore.Fields.zeros(surface_space) .+ FT(280.0) @test ClimaLand.surface_emissivity(model, Y, p) == emissivity - @test ClimaLand.surface_evaporative_scaling(model, Y, p) == FT(1) @test ClimaLand.surface_height(model, Y, p) == z_sfc @test ClimaLand.surface_albedo(model, Y, p) == PAR_albedo / 2 + NIR_albedo / 2 @test ClimaLand.surface_temperature(model, Y, p, t) == T_sfc - thermo_params = - LP.thermodynamic_parameters(model.parameters.earth_param_set) - ts_in = - Thermodynamics.PhaseEquil_pTq.( - thermo_params, - p.drivers.P, - p.drivers.T, - p.drivers.q, - ) - ρ_sfc = compute_ρ_sfc.(thermo_params, ts_in, T_sfc) - @test ClimaLand.surface_air_density( - model.boundary_conditions.top.atmos, - model, - Y, - p, - t, - T_sfc, - ) == ρ_sfc - - q_sat = - Thermodynamics.q_vap_saturation_generic.( - Ref(thermo_params), - T_sfc, - ρ_sfc, - Ref(Thermodynamics.Liquid()), - ) - g = LP.grav(model.parameters.earth_param_set) - M_w = LP.molar_mass_water(model.parameters.earth_param_set) - R = LP.gas_constant(model.parameters.earth_param_set) - ψ_sfc = p.soil.sfc_scratch - ClimaLand.Domains.linear_interpolation_to_surface!( - ψ_sfc, - p.soil.ψ, - coords.subsurface.z, - Δz_top, - ) - q_sfc = @. (q_sat * exp(g * ψ_sfc * M_w / (R * T_sfc))) - @test ClimaLand.surface_specific_humidity( - model, - Y, - p, - T_sfc, - ρ_sfc, - ) == q_sfc - conditions = ClimaLand.turbulent_fluxes( model.boundary_conditions.top.atmos, model, @@ -275,34 +226,35 @@ for FT in (Float32, Float64) ) computed_water_flux = p.soil.top_bc.water computed_energy_flux = p.soil.top_bc.heat - (; ν, θ_r, d_ds) = model.parameters - _D_vapor = FT(LP.D_vapor(model.parameters.earth_param_set)) - S_l_sfc = p.soil.sfc_scratch - ClimaLand.Domains.linear_interpolation_to_surface!( - S_l_sfc, - Soil.effective_saturation.(ν, Y.soil.ϑ_l, θ_r), - coords.subsurface.z, - Δz_top, - ) - τ_a = ClimaLand.Domains.top_center_to_surface( - @. (ν - p.soil.θ_l - Y.soil.θ_i)^(FT(5 / 2)) / ν - ) - f_ice = ClimaLand.Domains.top_center_to_surface( - @. effective_saturation(ν, Y.soil.θ_i, θ_r) / ( - effective_saturation(ν, Y.soil.θ_i, θ_r) + - effective_saturation(ν, Y.soil.ϑ_l, θ_r) - ) - ) - dsl = Soil.dry_soil_layer_thickness.(S_l_sfc, S_c, d_ds) - r_soil = @. dsl / (_D_vapor * τ_a) # [s\m] - r_ae = conditions.r_ae - expected_water_flux = @. FT(precip(t)) .+ - conditions.vapor_flux * (1 - f_ice) * r_ae / (r_soil + r_ae) + + expected_water_flux = @. FT(precip(t)) .+ conditions.vapor_flux_liq @test computed_water_flux == expected_water_flux - expected_energy_flux = @. R_n + - conditions.lhf * r_ae / (r_soil + r_ae) + - conditions.shf + expected_energy_flux = @. R_n + conditions.lhf + conditions.shf @test computed_energy_flux == expected_energy_flux + + # Test soil resistances for liquid water + # ϑ_sfc = range(θ_r-eps(FT), ν, 5) + # S_sfc = @. ClimaLand.Soil.effective_saturation(ν, ϑ_sfc, θ_r) + # S_i_sfc = range(FT(0), ν, 5) ./ ν; + # log10_K_sfc = @. log10(hydraulic_conductivity( hcm, K_sat, S_sfc) .+ eps(FT)) + # thermo_params = LP.thermodynamic_parameters(earth_param_set) + # ts_in = Thermodynamics.PhaseEquil_ρTq(thermo_params,FT(1.235), FT(285),FT(0.005)) + # fluxes = soil_turbulent_fluxes_at_a_point( + # parent(T_sfc)[1], + # parent(h_sfc)[1], + # parent(d_sfc)[1], + # S_i_sfc, + # hcm, + # K_sat, + # log10_K_sfc, + # ts_in, + # FT(3), + # atmos.h, + # FT(0), + # z_0m, + # z_0b, + # earth_param_set, + #) end end end diff --git a/test/standalone/Soil/soil_parameterizations.jl b/test/standalone/Soil/soil_parameterizations.jl index 277febf370..096b8ee156 100644 --- a/test/standalone/Soil/soil_parameterizations.jl +++ b/test/standalone/Soil/soil_parameterizations.jl @@ -184,26 +184,6 @@ for FT in (Float32, Float64) @test viscosity_factor.(T, parameters.γ, parameters.γT_ref) ≈ exp.(parameters.γ .* (T .- parameters.γT_ref)) - x = [0.0, 0.2, 0.4] - @test soil_tortuosity.(x, 0.0, 0.3) ≈ - [0.3^1.5, 0.1^2.5 / 0.3, eps(FT)^2.5 / 0.3] - x = [0.35, 0.25, 0.0, 0.1] - y = [0.0, 0.0, 0.1, 0.1 * 0.15 / 0.25] - @test dry_soil_layer_thickness.(x, 0.25, 0.1) ≈ y - @test soil_resistance( - θ_r + eps(FT), - θ_r, - parameters.hydrology_cm, - parameters.ν, - parameters.θ_r, - parameters.d_ds, - param_set, - ) ≈ - dry_soil_layer_thickness( - Soil.effective_saturation(ν, 2 * θ_r + eps(FT), θ_r), - hcm.S_c, - parameters.d_ds, - ) / FT(LP.D_vapor(param_set)) / soil_tortuosity(θ_r + eps(FT), θ_r, ν) end @testset "Brooks and Corey closure, FT = $FT" begin diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl index 28d12d6581..daf4034d80 100644 --- a/test/standalone/Vegetation/canopy_model.jl +++ b/test/standalone/Vegetation/canopy_model.jl @@ -236,13 +236,7 @@ import ClimaParams # check that this is updated correctly: # @test p.canopy.autotrophic_respiration.Ra == exp_tendency!(dY, Y, p, t0) - turb_fluxes = ClimaLand.Canopy.canopy_turbulent_fluxes( - canopy.atmos, - canopy, - Y, - p, - t0, - ) + turb_fluxes = ClimaLand.turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux @test p.canopy.energy.shf == turb_fluxes.shf @@ -731,13 +725,7 @@ end exp_tendency! = make_exp_tendency(canopy) dY = similar(Y) exp_tendency!(dY, Y, p, t0) - turb_fluxes = ClimaLand.Canopy.canopy_turbulent_fluxes( - canopy.atmos, - canopy, - Y, - p, - t0, - ) + turb_fluxes = ClimaLand.turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux @test p.canopy.energy.lhf == turb_fluxes.lhf @test p.canopy.energy.shf == turb_fluxes.shf