Skip to content

Commit

Permalink
photosynthesis_jl passed test
Browse files Browse the repository at this point in the history
  • Loading branch information
kongdd committed Oct 13, 2024
1 parent 9e9ef75 commit f0170fa
Show file tree
Hide file tree
Showing 11 changed files with 90 additions and 79 deletions.
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -29,3 +29,5 @@ docs/reference
__pycache__
data
examples

deps
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "BEPS"
uuid = "28258c4d-138d-4793-867a-264e32c782b1"
authors = ["Dongdong Kong <[email protected]>"]
version = "0.1.3"
version = "0.1.5"

[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Expand Down
11 changes: 10 additions & 1 deletion src/photosynthesis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ since transfer is going on only one side of a leaf.
(2PPFD < 1) && (PPFD = 0.0)
T_leaf_K = T_leaf + 273.13

λ = LAMBDA(T_leaf_p) # [J kg-1], ~2.5MJ kg-1
λ = leaf_lambda(T_leaf_p) # [J kg-1], ~2.5MJ kg-1
rᵥ = 1.0 / gb_w

ρₐ = cal_rho_a(T_leaf, ea) # [kg m-3]
Expand Down Expand Up @@ -133,6 +133,15 @@ end
exp(tprime * eact / (tref * rugc * t_lk))
end

function leaf_lambda(TK::T)::T where {T<:Real}
y = 3149_000.0 - 2370.0 * TK # J kg-1
# add heat of fusion for melting ice
if TK < 273.0
y += 333.0 # TODO: unit error, `y += 333_000.0`
end
return y
end

"""
If `Wj` or `Wc` are less than Rd then A would probably be less than 0. This
would yield a negative stomatal conductance.
Expand Down
9 changes: 0 additions & 9 deletions src/photosynthesis_helper.jl
Original file line number Diff line number Diff line change
Expand Up @@ -107,15 +107,6 @@ end
return rh
end

function LAMBDA(TK::T)::T where {T<:Real}
y = 3149000.0 - 2370.0 * TK # J kg-1
# add heat of fusion for melting ice
if TK < 273.0
y += 333.0 # TODO: unit error, `y += 333_000.0`
end
return y
end

# Function to calculate saturation vapor pressure function in mb
@fastmath function ES(t::T)::T where {T<:Real}
if t > 0.0
Expand Down
42 changes: 0 additions & 42 deletions test/debug-photo.qmd

This file was deleted.

4 changes: 4 additions & 0 deletions test/modules/modules.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
include("test-snowpack.jl")
include("test-photosynthesis.jl")
include("test-sensible_heat.jl")
include("test-aerodynamic_conductance.jl")
23 changes: 23 additions & 0 deletions test/modules/test-aerodynamic_conductance.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
@testset "aerodynamic_conductance" begin
canopyh_o = 2.0
canopyh_u = 0.2
height_wind_sp = 2.0
clumping = 0.8
Ta = 20.0
wind_sp = 2.0
GH_o = 100.0
pai_o = 4.0
pai_u = 2.0

ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
pai_o, pai_u)
r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)

ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
clang.aerodynamic_conductance_c(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
pai_o, pai_u)
r2 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)

@test maximum(abs.(values(r1) .- values(r2))) <= 1e-8
end
30 changes: 30 additions & 0 deletions test/modules/test-photosynthesis.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
@testset "photosynthesis_jl" begin
ea = 0.9879918003910609 # kPa
gb_w = 0.60 # m s-1
VCmax25 = 124.13532190023277 # μmol m-2 s-1

β_soil = 0.9514687124755321
g0_w = 0.0175
g1_w = 8.0

cii = 266.0
LH_leaf = 68.04830184404008

## 白天
Tleaf = 16.0
Tleaf_c = 16.5
Rs_leaf = 378.15760935295305

r1 = clang.photosynthesis_c(Tleaf, Rs_leaf, ea, gb_w, VCmax25, β_soil, g0_w, g1_w, cii, Tleaf_c, LH_leaf)
r2 = photosynthesis_jl(Tleaf, Rs_leaf, ea, gb_w, VCmax25, β_soil, g0_w, g1_w, cii, Tleaf_c, LH_leaf)
@test all(isapprox.(r1, r2, rtol=1e-7))

## 夜间
Tleaf = 16.0
Tleaf_c = 16.5
Rs_leaf = 0.0

r1 = clang.photosynthesis_c(Tleaf, Rs_leaf, ea, gb_w, VCmax25, β_soil, g0_w, g1_w, cii, Tleaf_c, LH_leaf)
r2 = photosynthesis_jl(Tleaf, Rs_leaf, ea, gb_w, VCmax25, β_soil, g0_w, g1_w, cii, Tleaf_c, LH_leaf)
@test all(isapprox.(r1, r2, rtol=1e-8))
end
17 changes: 17 additions & 0 deletions test/modules/test-sensible_heat.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
@testset "sensible_heat_jl" begin
T_leaf = Leaf(20.0)
T_ground = 15.0
T_air = 22.0
rh_air = 80.0

Gheat = Leaf(100.0)
Gheat_g = 120.0
lai = Leaf(2.0)

r1 = clang.sensible_heat_c(T_leaf, T_ground, T_air, rh_air,
Gheat, Gheat_g, lai)

r2 = sensible_heat_jl(T_leaf, T_ground, T_air, rh_air,
Gheat, Gheat_g, lai)
all(r1 .≈ r2)
end
5 changes: 3 additions & 2 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@ using Test
using BEPS
using BEPS: path_proj

include("modules/test-snowpack.jl")
# include("test-beps_main.jl")
include("test-clang.jl")
include("modules/modules.jl")

include("test-soil.jl")
include("test-clang.jl")
24 changes: 0 additions & 24 deletions test/test-clang.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,27 +20,3 @@ end
@test maximum(abs.(coef1 .- coef2)) == 0
end
end

@testset "aerodynamic_conductance" begin
canopyh_o = 2.0
canopyh_u = 0.2
height_wind_sp = 2.0
clumping = 0.8
Ta = 20.0
wind_sp = 2.0
GH_o = 100.0
pai_o = 4.0
pai_u = 2.0

ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
pai_o, pai_u)
r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)

ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
clang.aerodynamic_conductance_c(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
pai_o, pai_u)
r2 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)

@test maximum(abs.(values(r1) .- values(r2))) <= 1e-8
end

0 comments on commit f0170fa

Please sign in to comment.