Skip to content

Commit

Permalink
rename symbols
Browse files Browse the repository at this point in the history
  • Loading branch information
kongdd committed Oct 13, 2024
1 parent 489fdd3 commit ec241e8
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 48 deletions.
15 changes: 7 additions & 8 deletions src/inter_prg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ The inter-module function between main program and modules
"""
function inter_prg_jl(
jday::Int, hour::Int,
lai::T, clumping::T, param::Vector{T}, meteo::Met, CosZs::T,
lai::T, Ω::T, param::Vector{T}, meteo::Met, CosZs::T,
var_o::Vector{T}, var_n::Vector{T}, soil::AbstractSoil,
Ra::Radiation,
mid_res::Results, mid_ET::OutputET, var::InterTempVars; kw...) where {T}
Expand Down Expand Up @@ -49,7 +49,7 @@ function inter_prg_jl(
Tcu = 0.0

# /***** Vcmax-Nitrogen calculations,by G.Mo,Apr. 2011 *****/
Vcmax_sunlit, Vcmax_shaded = VCmax(lai, clumping, CosZs, param)
Vcmax_sunlit, Vcmax_shaded = VCmax(lai, Ω, CosZs, param)

# /***** LAI calculation module, by B. Chen *****/
lai_o = lai < 0.1 ? 0.1 : lai
Expand All @@ -66,7 +66,7 @@ function inter_prg_jl(
stem_u = param[9+1] * 0.2 # parameter[9].LAI max understory

# lai2: separate lai into sunlit and shaded portions
lai2(clumping, CosZs, stem_o, stem_u, lai_o, lai_u, LAI, PAI)
lai2(Ω, CosZs, stem_o, stem_u, lai_o, lai_u, LAI, PAI)

# /***** Initialization of this time step *****/
Rs = meteo.Srad
Expand Down Expand Up @@ -136,14 +136,13 @@ function inter_prg_jl(

# set!(m_snow_pre, 0.0);
depth_snow = snowpack_stage1_jl(Ta, prcp,
lai_o, lai_u,
clumping,
lai_o, lai_u, Ω,
m_snow_pre, m_snow, f_snow, A_snow,
depth_snow, ρ_snow,
α_v_sw, α_n_sw) # by X. Luo

# /***** Rain fall stage 1 by X. Luo *****/
var.r_rain_g[k] = rainfall_stage1_jl(Ta, prcp, f_water, m_water, m_water_pre, lai_o, lai_u, clumping)
var.r_rain_g[k] = rainfall_stage1_jl(Ta, prcp, f_water, m_water, m_water_pre, lai_o, lai_u, Ω)

if (soil.θ_prev[2] < soil.θ_vwp[2] * 0.5)
α_g = α_dry
Expand Down Expand Up @@ -173,7 +172,7 @@ function inter_prg_jl(
num = num + 1
# /***** Aerodynamic_conductance module by G.Mo *****/
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, GH_o,
aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, Ω, Ta, wind, GH_o,
lai_o + stem_o, lai_u + stem_u)

init_leaf_dbl2(Gh,
Expand All @@ -192,7 +191,7 @@ function inter_prg_jl(
# /***** Net radiation at canopy and leaf level module by X.Luo *****/
radiation_o, radiation_u, radiation_g = netRadiation_jl(Rs, CosZs, Tco, Tcu, temp_grd,
lai_o, lai_u, lai_o + stem_o, lai_u + stem_u, PAI,
clumping, Ta, RH,
Ω, Ta, RH,
α_v_sw[], α_n_sw[],
percArea_snow_o, percArea_snow_u,
f_snow.g,
Expand Down
33 changes: 17 additions & 16 deletions src/snowpack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@ end
# Arguments
- `m_snow`: m_snow
- `mw`: m_water
- `depth_snow`: snow depth
- `depth_water`: water depth
- `z_snow`: snow depth
- `z_water`: water depth
*reference variables*
- m_snow_pre
Expand All @@ -44,7 +44,7 @@ function snowpack_stage1_jl(Tair::Float64, prcp::Float64,
m_snow::Layer3{Float64},
perc_snow::Layer3{Float64},
area_snow::Layer2{Float64},
depth_snow::Float64,
z_snow::Float64,
ρ_snow::Ref{Float64},
albedo_v_snow::Ref{Float64}, albedo_n_snow::Ref{Float64})

Expand Down Expand Up @@ -92,12 +92,13 @@ function snowpack_stage1_jl(Tair::Float64, prcp::Float64,
m_snow.g = max(0.0, m_snow_pre.g + δ_zs * ρ_new) # [kg m-2]

if δ_zs > 0
ρ_snow[] = (ρ_snow[] * depth_snow + ρ_new * δ_zs) / (depth_snow + δ_zs) # 计算混合密度
# 可能导致了ρ_snow的极低。
ρ_snow[] = (ρ_snow[] * z_snow + ρ_new * δ_zs) / (z_snow + δ_zs) # 计算混合密度
else
ρ_snow[] = (ρ_snow[] - 250) * exp(-0.001 * kstep / 3600.0) + 250.0
end

depth_snow = m_snow.g > 0 ? m_snow.g / ρ_snow[] : 0.0
z_snow = m_snow.g > 0 ? m_snow.g / ρ_snow[] : 0.0
perc_snow.g = min(m_snow.g / (0.05 * ρ_snow[]), 1.0) # [m],认为雪深50cm时,perc_snow=1

if snowrate_o > 0
Expand All @@ -107,7 +108,7 @@ function snowpack_stage1_jl(Tair::Float64, prcp::Float64,
albedo_v_snow[] = albedo_v_new
albedo_n_snow[] = albedo_n_new
end
min(depth_snow, 10.0) # 雪深过高,限制为10m即可
min(z_snow, 10.0) # 雪深过高,限制为10m即可
end


Expand All @@ -121,11 +122,11 @@ end
# 热量释放用于融雪(Tsnow -> 0)
# - pos: 融化
# - neg: 冻结
function cal_melt(depth_snow::T, ρ_snow::T, Tsnow::T) where {T<:Real}
function cal_melt(z_snow::T, ρ_snow::T, Tsnow::T) where {T<:Real}
dT = Tsnow - 0
cp_ice = 2228.261 # J Kg-1 K-1
λ_fusion = 3.34 * 1000000 # J Kg-1
m = depth_snow * ρ_snow # [kg m-2]
m = z_snow * ρ_snow # [kg m-2]
E = m * dT * cp_ice # 当前的雪融化,需要这么多能量,
return E / λ_fusion # m_ice, -> kg
end
Expand All @@ -143,14 +144,14 @@ It is assumed sublimation happens before the melting and freezing process.
- depth: [m]
"""
function snowpack_stage3_jl(Tair::Float64, Tsnow::Float64, Tsnow_last::Float64, ρ_snow::Float64,
depth_snow::Float64, depth_water::Float64, m_snow::Layer3{Float64})
z_snow::Float64, z_water::Float64, m_snow::Layer3{Float64})

zs_sup = depth_snow # already considered sublimation
zs_sup = z_snow # already considered sublimation
ms_sup = m_snow.g

Δm = cal_melt(zs_sup, ρ_snow, Tsnow) # [kg m-2]
con_melt = Tsnow > 0 && Tsnow_last <= 0 && ms_sup > 0
con_frozen = Tsnow <= 0 && Tsnow_last > 0 && depth_water > 0
con_frozen = Tsnow <= 0 && Tsnow_last > 0 && z_water > 0

ms_melt = 0.0
mw_frozen = 0.0
Expand All @@ -163,16 +164,16 @@ function snowpack_stage3_jl(Tair::Float64, Tsnow::Float64, Tsnow_last::Float64,
elseif 0.02 < zs_sup <= 0.05
# case 2 depth of snow > 0.02 < 0.05 m
con_melt && (ms_melt = min(Δm, ms_sup))
con_frozen && (mw_frozen = min(-Δm, depth_water * ρ_w))
con_frozen && (mw_frozen = min(-Δm, z_water * ρ_w))
elseif zs_sup > 0.05
# _z = max(zs_sup*0.02, 0.02) # TODO: fix释放热量的深度
# _Δm = cal_melt(_z, Tsnow)
con_melt && (ms_melt = min(0.02Δm, ms_sup))
con_frozen && (mw_frozen = min(-0.02Δm, depth_water * ρ_w))
con_frozen && (mw_frozen = min(-0.02Δm, z_water * ρ_w))
end

m_snow.g = max(0.0, m_snow.g - ms_melt + mw_frozen) # ground snow
depth_snow = max(0.0, zs_sup + (mw_frozen - ms_melt) / ρ_snow)
depth_water = max(0.0, depth_water + (ms_melt - mw_frozen) / ρ_w)
depth_snow, depth_water
z_snow = max(0.0, zs_sup + (mw_frozen - ms_melt) / ρ_snow)
z_water = max(0.0, z_water + (ms_melt - mw_frozen) / ρ_w)
z_snow, z_water
end
48 changes: 24 additions & 24 deletions src/surface_temperature.jl
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
"""
- cp: capacity_heat, [J kg-1 K-1]
"""
function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_water::FT,
function surface_temperature_jl(T_air::FT, rh_air::FT, z_snow::FT, z_water::FT,
cp_soil1::FT, cp_soil0::FT, Gheat_g::FT,
depth_soil1::FT, ρ_snow::FT, tempL_u::FT, Rn_g::FT,
z_soil1::FT, ρ_snow::FT, tempL_u::FT, Rn_g::FT,
E_soil::FT, E_water_g::FT, E_snow_g::FT, λ_soil1::FT,
perc_snow_g::FT, G_soil1::FT,
T_ground_last::FT,
Expand Down Expand Up @@ -32,10 +32,10 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat
T_snow2::FT = 0.0
G::FT = 0.0

if depth_snow <= 0.02
if z_snow <= 0.02
ttt = cp_soil1 * 0.02 / length_step
T_ground = (T_ground_last * ttt * ra_g * depth_soil1 + Gg * ra_g * depth_soil1 + ρₐ * cp * T_air * depth_soil1 + ra_g * λ_soil1 * T_soil1_last) /
(ρₐ * cp * depth_soil1 + ra_g * λ_soil1 + ttt * ra_g * depth_soil1)
T_ground = (T_ground_last * ttt * ra_g * z_soil1 + Gg * ra_g * z_soil1 + ρₐ * cp * T_air * z_soil1 + ra_g * λ_soil1 * T_soil1_last) /
(ρₐ * cp * z_soil1 + ra_g * λ_soil1 + ttt * ra_g * z_soil1)
T_ground = clamp(T_ground, T_ground_last - 25, T_ground_last + 25)

T_any0 = T_ground
Expand All @@ -44,39 +44,39 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat
T_snow1 = T_any0
T_snow2 = T_any0

G = 2 * λ_soil1 * (T_any0 - T_soil1_last) / depth_soil1
G = 2 * λ_soil1 * (T_any0 - T_soil1_last) / z_soil1
G = clamp(G, -100.0, 100.0)

elseif depth_snow > 0.02 && depth_snow <= 0.05
elseif z_snow > 0.02 && z_snow <= 0.05
ttt = cp_soil1 * 0.02 / length_step # for soil fraction part

T_soil0 = (T_soil0_last * ttt * ra_g * depth_soil1 + Gg * ra_g * depth_soil1 + ρₐ * cp * T_air * depth_soil1 + 2 * ra_g * λ_soil1 * T_soil1_last) /
(ρₐ * cp * depth_soil1 + 2 * ra_g * λ_soil1 + ttt * ra_g * depth_soil1)
T_soil0 = (T_soil0_last * ttt * ra_g * z_soil1 + Gg * ra_g * z_soil1 + ρₐ * cp * T_air * z_soil1 + 2 * ra_g * λ_soil1 * T_soil1_last) /
(ρₐ * cp * z_soil1 + 2 * ra_g * λ_soil1 + ttt * ra_g * z_soil1)

T_soil0 = clamp(T_soil0, T_air - 25.0, T_air + 25.0)

ttt = cp_ice * ρ_snow * depth_snow / length_step # for snow part
T_snow = (T_snow_last * ttt * ra_g * depth_snow + Gg * ra_g * depth_snow + ρₐ * cp * tempL_u * depth_snow + ra_g * κ_snow * T_any0_last) /
(ρₐ * cp * depth_snow + ra_g * κ_snow + ttt * ra_g * depth_snow)
ttt = cp_ice * ρ_snow * z_snow / length_step # for snow part
T_snow = (T_snow_last * ttt * ra_g * z_snow + Gg * ra_g * z_snow + ρₐ * cp * tempL_u * z_snow + ra_g * κ_snow * T_any0_last) /
(ρₐ * cp * z_snow + ra_g * κ_snow + ttt * ra_g * z_snow)

T_snow = clamp(T_snow, T_air - 25.0, T_air + 25.0)

ttt = (λ_soil1 * T_soil1_last / depth_soil1 + T_snow * κ_snow + 0.02 * cp_soil1 / length_step * T_any0_last) /
(λ_soil1 / depth_soil1 + κ_snow / depth_snow + 0.02 * cp_soil1 / length_step)
ttt = (λ_soil1 * T_soil1_last / z_soil1 + T_snow * κ_snow + 0.02 * cp_soil1 / length_step * T_any0_last) /
(λ_soil1 / z_soil1 + κ_snow / z_snow + 0.02 * cp_soil1 / length_step)
T_any0 = T_soil0 * (1 - perc_snow_g) + ttt * perc_snow_g

G_snow = κ_snow / (depth_snow + 0.5 * depth_soil1) * (T_snow - T_soil1_last)
G_soil = G_snow * (T_any0 - T_soil1_last) / depth_soil1
G_snow = κ_snow / (z_snow + 0.5 * z_soil1) * (T_snow - T_soil1_last)
G_soil = G_snow * (T_any0 - T_soil1_last) / z_soil1

G = G_snow * perc_snow_g + G_soil * (1 - perc_snow_g)
G = clamp(G, -100.0, 100.0)

# starting to melt
if T_snow > 0.0 && T_snow_last <= 0.0 && depth_snow > 0.0
if T_snow > 0.0 && T_snow_last <= 0.0 && z_snow > 0.0
T_snow = 0.0
end
# starting to frozen
if T_snow < 0.0 && T_snow_last >= 0.0 && depth_water > 0.0
if T_snow < 0.0 && T_snow_last >= 0.0 && z_water > 0.0
T_snow = 0.0
end

Expand All @@ -86,7 +86,7 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat

T_snow1 = T_snow
T_snow2 = T_snow
elseif depth_snow > 0.05
elseif z_snow > 0.05
ttt = cp_ice * ρ_snow * 0.02 / length_step

T_snow = (T_snow_last * ttt * ra_g * 0.04 + Gg * ra_g * 0.02 + ρₐ * cp * T_air * 0.04 + ra_g * κ_snow * T_snow1_last) /
Expand All @@ -96,18 +96,18 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat
G_snow = κ_snow * (T_snow - T_snow1_last) / 0.04
G = clamp(G_snow, -100.0, 100.0)

G_snow1 = κ_snow * (T_snow1_last - T_snow2_last) / (depth_snow - 0.02)
G_snow1 = κ_snow * (T_snow1_last - T_snow2_last) / (z_snow - 0.02)
T_snow1 = T_snow1_last + ((G - G_snow1) / (cp_ice * ρ_snow * 0.02)) * length_step
G_snow2 = (T_snow2_last - T_any0_last) / ((0.5 * (depth_snow - 0.04) / κ_snow) + (0.02 / λ_soil1))
T_snow2 = T_snow2_last + ((G_snow1 - G_snow2) / (cp_ice * ρ_snow * (depth_snow - 0.04))) * length_step
G_snow2 = (T_snow2_last - T_any0_last) / ((0.5 * (z_snow - 0.04) / κ_snow) + (0.02 / λ_soil1))
T_snow2 = T_snow2_last + ((G_snow1 - G_snow2) / (cp_ice * ρ_snow * (z_snow - 0.04))) * length_step

T_any0 = T_any0_last + ((G_snow2 - G_soil1) / (cp_soil0 * 0.02)) * length_step
T_soil0 = T_any0

if T_snow > 0.0 && T_snow_last <= 0.0 && depth_snow > 0.0
if T_snow > 0.0 && T_snow_last <= 0.0 && z_snow > 0.0
T_snow = 0
end
if T_snow < 0.0 && T_snow_last >= 0.0 && depth_water > 0.0
if T_snow < 0.0 && T_snow_last >= 0.0 && z_water > 0.0
T_snow = 0
end
T_ground = T_snow
Expand Down

0 comments on commit ec241e8

Please sign in to comment.