Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

extraction fix #758

Merged
merged 1 commit into from
Sep 12, 2024
Merged

extraction fix #758

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/src/APIs/canopy/PlantHydraulics.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ ClimaLand.PlantHydraulics.augmented_liquid_fraction
ClimaLand.PlantHydraulics.water_retention_curve
ClimaLand.PlantHydraulics.inverse_water_retention_curve
ClimaLand.PlantHydraulics.root_water_flux_per_ground_area!
ClimaLand.PlantHydraulics.flux
ClimaLand.PlantHydraulics.water_flux
ClimaLand.PlantHydraulics.hydraulic_conductivity
```

Expand Down
17 changes: 14 additions & 3 deletions experiments/long_runs/land.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ import ClimaAnalysis
import ClimaAnalysis.Visualize as viz
import ClimaUtilities

import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
import ClimaUtilities.TimeVaryingInputs:
TimeVaryingInput, LinearInterpolation, PeriodicCalendar
import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
import ClimaUtilities.Regridders: InterpolationsRegridder
import ClimaUtilities.ClimaArtifacts: @clima_artifact
Expand All @@ -47,7 +48,7 @@ using Dates
import NCDatasets

const FT = Float64;

time_interpolation_method = LinearInterpolation(PeriodicCalendar())
regridder_type = :InterpolationsRegridder
context = ClimaComms.context()
device = ClimaComms.device()
Expand Down Expand Up @@ -83,6 +84,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
method = time_interpolation_method,
)

snow_precip = TimeVaryingInput(
Expand All @@ -92,6 +94,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
method = time_interpolation_method,
)

u_atmos = TimeVaryingInput(
Expand All @@ -100,20 +103,23 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)
q_atmos = TimeVaryingInput(
joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
"q",
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)
P_atmos = TimeVaryingInput(
joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
"sp",
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)

T_atmos = TimeVaryingInput(
Expand All @@ -122,6 +128,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)
h_atmos = FT(10)

Expand All @@ -145,6 +152,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
method = time_interpolation_method,
)
LW_d = TimeVaryingInput(
joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
Expand All @@ -153,6 +161,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
method = time_interpolation_method,
)

function zenith_angle(
Expand Down Expand Up @@ -489,6 +498,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
file_reader_kwargs = (;
preprocess_func = (data) -> data > 0.05 ? data : 0.0,
),
method = time_interpolation_method,
)
ai_parameterization =
Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
Expand Down Expand Up @@ -613,6 +623,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
start_date;
output_writer = nc_writer,
output_vars = :short,
average_period = :monthly,
)

diagnostic_handler =
Expand All @@ -627,7 +638,7 @@ end
function setup_and_solve_problem(; greet = false)

t0 = 0.0
tf = 60 * 60.0 * 24 * 14
tf = 60 * 60.0 * 24 * 365
Δt = 900.0
nelements = (101, 15)
if greet
Expand Down
15 changes: 12 additions & 3 deletions experiments/long_runs/soil.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ using ClimaAnalysis
import ClimaAnalysis.Visualize as viz
using ClimaUtilities

import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
import ClimaUtilities.TimeVaryingInputs:
TimeVaryingInput, LinearInterpolation, PeriodicCalendar
import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
import ClimaUtilities.Regridders: InterpolationsRegridder
import ClimaUtilities.ClimaArtifacts: @clima_artifact
Expand All @@ -45,7 +46,7 @@ using Dates
import NCDatasets

const FT = Float64;

time_interpolation_method = LinearInterpolation(PeriodicCalendar())
regridder_type = :InterpolationsRegridder
context = ClimaComms.context()
device = ClimaComms.device()
Expand Down Expand Up @@ -81,6 +82,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
method = time_interpolation_method,
)

snow_precip = TimeVaryingInput(
Expand All @@ -90,6 +92,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
method = time_interpolation_method,
)

u_atmos = TimeVaryingInput(
Expand All @@ -98,20 +101,23 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)
q_atmos = TimeVaryingInput(
joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
"q",
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)
P_atmos = TimeVaryingInput(
joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
"sp",
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)

T_atmos = TimeVaryingInput(
Expand All @@ -120,6 +126,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
surface_space;
reference_date = start_date,
regridder_type,
method = time_interpolation_method,
)
h_atmos = FT(10)

Expand All @@ -143,6 +150,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
method = time_interpolation_method,
)
LW_d = TimeVaryingInput(
joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
Expand All @@ -151,6 +159,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
reference_date = start_date,
regridder_type,
file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
method = time_interpolation_method,
)

function zenith_angle(
Expand Down Expand Up @@ -430,7 +439,7 @@ end
function setup_and_solve_problem(; greet = false)

t0 = 0.0
tf = 60 * 60.0 * 24 * 340
tf = 60 * 60.0 * 24 * 365
Δt = 900.0
nelements = (101, 15)
if greet
Expand Down
14 changes: 9 additions & 5 deletions src/integrated/soil_canopy_model.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -263,18 +263,22 @@ function make_update_boundary_fluxes(
area_index,
land.canopy.hydraulics.compartment_labels[1],
)
# Note that we model the flux between each soil layer and the canopy as:
# Flux = -K_eff x [(ψ_canopy - ψ_soil)/(z_canopy - z_soil) + 1], where
# K_eff = K_soil K_canopy /(K_canopy + K_soil)

# Note that in `PrescribedSoil` mode, we compute the flux using K_soil = K_plant(ψ_soil)
# and K_canopy = K_plant(ψ_canopy). In `PrognosticSoil` mode here, we compute the flux using
# K_soil = K_soil(ψ_soil) and K_canopy = K_plant(ψ_canopy).

@. p.root_extraction =
above_ground_area_index *
PlantHydraulics.flux(
PlantHydraulics.water_flux(
z,
land.canopy.hydraulics.compartment_midpoints[1],
p.soil.ψ,
p.canopy.hydraulics.ψ.:1,
PlantHydraulics.hydraulic_conductivity(
conductivity_model,
p.soil.ψ,
),
p.soil.K,
PlantHydraulics.hydraulic_conductivity(
conductivity_model,
p.canopy.hydraulics.ψ.:1,
Expand Down
2 changes: 1 addition & 1 deletion src/standalone/Vegetation/Canopy.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -524,7 +524,7 @@ function ClimaLand.make_update_aux(
# Compute the flux*area between the current compartment `i`
# and the compartment above.
@. fa.:($$i) =
PlantHydraulics.flux(
PlantHydraulics.water_flux(
hydraulics.compartment_midpoints[i],
hydraulics.compartment_midpoints[ip1],
ψ.:($$i),
Expand Down
47 changes: 31 additions & 16 deletions src/standalone/Vegetation/PlantHydraulics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ import ClimaLand:
name
export PlantHydraulicsModel,
AbstractPlantHydraulicsModel,
flux,
water_flux,
effective_saturation,
augmented_liquid_fraction,
water_retention_curve,
Expand Down Expand Up @@ -328,7 +328,7 @@ end


"""
flux(
water_flux(
z1,
z2,
ψ1,
Expand All @@ -337,18 +337,25 @@ end
K2,
) where {FT}

Computes the water flux given the absolute potential (pressure/(ρg))
at the center of the two compartments z1 and z2,
and the conductivity along
the flow path between these two points.
Computes the water flux given the absolute potential ψ (pressure/(ρg))
and the conductivity K (m/s) at the center of the two layers
with midpoints z1 and z2.

We currently assuming an arithmetic
mean for mean K_sat between the two points (Bonan, 2019; Zhu, 2008)
to take into account the change in K_sat halfway
between z1 and z2; this is incorrect for compartments of differing sizes.
We currently assuming a harmonic
mean for effective conducticity between the two layers
(see CLM Technical Documentation).

To account for different path lengths in the two compartments Δz1 and
Δz2, we would require the following conductance k (1/s)
k_eff = K1/Δz1*K2/Δz2/(K1/Δz1+K2/Δz2)
and a water flux of
F = -k_eff * (ψ1 +z1 - ψ2 - z2) (m/s).

This currently assumes the path lengths are equal.
"""
function flux(z1, z2, ψ1, ψ2, K1, K2)
flux = -(K1 + K2) / 2 * ((ψ2 - ψ1) / (z2 - z1) + 1)
function water_flux(z1::FT, z2::FT, ψ1::FT, ψ2::FT, K1::FT, K2::FT) where {FT}
K_eff = K1 * K2 / max(K1 + K2, eps(FT))
flux = -K_eff * ((ψ2 - ψ1) / (z2 - z1) + 1)
return flux # (m/s)
end

Expand Down Expand Up @@ -580,8 +587,16 @@ end
) where {FT}

A method which computes the water flux between the soil and the stem, via the roots,
and multiplied by the RAI, in the case of a model running without an integrated
soil model.
and multiplied by the RAI, in the case of a model running without a prognostic
soil model:

Flux = -K_eff x [(ψ_stem - ψ_soil)/(z_stem - z_soil) + 1], where
K_eff = K_soil K_stem /(K_stem + K_soil)

Note that in `PrescribedSoil` mode, we compute the flux using K_soil = K_plant(ψ_soil)
and K_stem = K_plant(ψ_stem). In `PrognosticSoil` mode, we compute the flux using
K_soil = K_soil(ψ_soil) and K_stem = K_plant(ψ_stem). The latter is a better model, but
our `PrescribedSoil` struct does not store K_soil, only ψ_soil.

The returned flux is per unit ground area. This assumes that the stem compartment
is the first element of `Y.canopy.hydraulics.ϑ_l`.
Expand Down Expand Up @@ -610,7 +625,7 @@ function root_water_flux_per_ground_area!(

if i != n_root_layers
@. fa +=
flux(
water_flux(
root_depths[i],
model.compartment_midpoints[1],
ψ_soil,
Expand All @@ -623,7 +638,7 @@ function root_water_flux_per_ground_area!(
above_ground_area_index
else
@. fa +=
flux(
water_flux(
root_depths[i],
model.compartment_midpoints[1],
ψ_soil,
Expand Down
Loading
Loading