patch

Author: kongdd

Diff for: src/BEPS/Base/VCmax.jl

@@ -1,4 +1,6 @@
-function VCmax(lai, clumping, CosZs, param)
+function VCmax(lai::FT, clumping::FT, CosZs::FT, param::AbstractVector{FT})
+  CosZs <= 0 && return 0.0, 0.0 # 光合仅发生在白天
+  
   # parameters for Vcmax-Nitrogen calculations
   G_theta = 0.5 # assuming a spherical leaf angle distribution
   K = G_theta * clumping / CosZs

Diff for: src/BEPS/evaporation_canopy.jl

@@ -1,6 +1,7 @@
 function evaporation_canopy_jl(T_leaf::Leaf, Ta::Float64, RH::Float64,
   Gwater::Leaf, lai::Leaf,
-  perc_water_o::Float64, perc_water_u::Float64, perc_snow_o::Float64, perc_snow_u::Float64)
+  perc_water_o::Float64, perc_water_u::Float64, 
+  perc_snow_o::Float64, perc_snow_u::Float64)
 
   # LHw = Leaf()  # latent heat from leaves W/m2, caused by evaporation of intercepted rain
   # LHs = Leaf()  # latent heat from leaves W/m2, caused by evaporation of intercepted snow

Diff for: src/BEPS/evaporation_soil.jl

@@ -8,7 +8,7 @@ function evaporation_soil_jl(Tair::FT, Tg::FT, RH::FT, Rn_g::FT, Gheat_g::FT,
   perc_snow_g::Ref{FT}, 
   depth_water, depth_snow,
   # depth_water::Ref{FT}, depth_snow::Ref{FT},
-  mass_water_g::FT, mass_snow_g::Ref{FT},
+  mass_water_g::FT, mass_snow::CanopyLayer{FT},
   ρ_snow::FT, swc_g::FT, porosity_g::FT) where {FT<:Real}
 
   met = meteo_pack_jl(Tg, RH)
@@ -18,7 +18,7 @@ function evaporation_soil_jl(Tair::FT, Tg::FT, RH::FT, Rn_g::FT, Gheat_g::FT,
   Gwater_g = 1.0 / (4.0 * exp(8.2 - 4.2 * swc_g / porosity_g))
   # kstep = kstep # 360s, 6min
 
-  perc_snow_g[] = depth_snow > 0.02 ? 1.0 : mass_snow_g[] / (0.025 * ρ_snow)
+  perc_snow_g[] = depth_snow > 0.02 ? 1.0 : mass_snow.g / (0.025 * ρ_snow)
   perc_snow_g[] = clamp(perc_snow_g[], 0.0, 1.0)
 
   if depth_water > 0.0 && depth_snow == 0.0
@@ -46,11 +46,11 @@ function evaporation_soil_jl(Tair::FT, Tg::FT, RH::FT, Rn_g::FT, Gheat_g::FT,
 
   Esoil_g = max(-0.002 / kstep, Esoil_g)
   if Esoil_g > 0.0
-    Esoil_g = min(Esoil_g, mass_snow_g[] / kstep)
+    Esoil_g = min(Esoil_g, mass_snow.g / kstep)
   end
 
-  mass_snow_g[] = max(mass_snow_g[] - Esoil_g * kstep, 0)
-  depth_snow = mass_snow_g[] > 0.0 ? depth_snow - (Esoil_g / ρ_snow) * kstep : 0.0
+  mass_snow.g = max(mass_snow.g - Esoil_g * kstep, 0)
+  depth_snow = mass_snow.g > 0.0 ? depth_snow - (Esoil_g / ρ_snow) * kstep : 0.0
 
   if depth_water > 0.0 || depth_snow > 0.0
     Esoil = 0.0

Diff for: src/BEPS/snowpack.jl

@@ -8,7 +8,8 @@
 """
 function snowpack_stage1_jl(temp_air::Float64, prcp::Float64,
   # ms_o_last::Float64, ms_u_last::Float64, ms_g_last::Float64, 
-  ms_o::Ref{Float64}, ms_u::Ref{Float64}, ms_g::Ref{Float64},
+  ms::CanopyLayer{Float64},
+  # ms_o::Ref{Float64}, ms_u::Ref{Float64}, ms_g::Ref{Float64},
   lai_o::Float64, lai_u::Float64, clumping::Float64, area_snow_o::Ref{Float64}, area_snow_u::Ref{Float64},
   # perc_snow_o::Ref{Float64}, perc_snow_u::Ref{Float64}, perc_snow_g::Ref{Float64},
   ρ_snow::Ref{Float64}, d_s, albedo_v_snow::Ref{Float64}, albedo_n_snow::Ref{Float64})
@@ -22,54 +23,54 @@ function snowpack_stage1_jl(temp_air::Float64, prcp::Float64,
   albedo_v_Newsnow = 0.94
   albedo_n_Newsnow = 0.8
 
-  ms_o_last = ms_o[]
-  ms_u_last = ms_u[]
-  ms_g_last = ms_g[]
+  ms_o_last = ms.o
+  ms_u_last = ms.u
+  ms_g_last = ms.g
 
   snowrate = temp_air > 0 ? 0 : prcp * ρ_w / ρ_new_snow
   snowrate_o = 0.0
 
   if temp_air < 0
     snowrate_o = snowrate
-    ms_o[] = ms_o_last + snowrate_o * kstep * ρ_new_snow * (1 - exp(-lai_o * clumping))
-    perc_snow_o = ms_o[] / massMax_snow_o
+    ms.o = ms_o_last + snowrate_o * kstep * ρ_new_snow * (1 - exp(-lai_o * clumping))
+    perc_snow_o = ms.o / massMax_snow_o
     perc_snow_o = clamp(perc_snow_o, 0, 1)
     area_snow_o[] = perc_snow_o * areaMax_snow_o
-    massStep_snow_o = ms_o[] - ms_o_last
+    massStep_snow_o = ms.o - ms_o_last
 
     snowrate_u = max(0, snowrate_o - (massStep_snow_o) / ρ_new_snow / kstep)
-    ms_u[] = ms_u_last + snowrate_u * kstep * ρ_new_snow * (1 - exp(-lai_u * clumping))
-    perc_snow_u = ms_u[] / massMax_snow_u
+    ms.u = ms_u_last + snowrate_u * kstep * ρ_new_snow * (1 - exp(-lai_u * clumping))
+    perc_snow_u = ms.u / massMax_snow_u
     perc_snow_u = clamp(perc_snow_u, 0, 1)
     area_snow_u[] = perc_snow_u * areaMax_snow_u
-    massStep_snow_u = ms_u[] - ms_u_last
+    massStep_snow_u = ms.u - ms_u_last
 
     snowrate_g = max(0, snowrate_u - (massStep_snow_u) / ρ_new_snow / kstep)
     δ_d_s = snowrate_g * kstep
   else
-    ms_o[] = ms_o_last
-    perc_snow_o = clamp(ms_o[] / massMax_snow_o, 0, 1)
+    ms.o = ms_o_last
+    perc_snow_o = clamp(ms.o / massMax_snow_o, 0, 1)
     area_snow_o[] = area_snow_o[]
 
-    ms_u[] = ms_u_last
-    perc_snow_u = ms_u[] / massMax_snow_u
+    ms.u = ms_u_last
+    perc_snow_u = ms.u / massMax_snow_u
     perc_snow_u = clamp(perc_snow_u, 0, 1)
     area_snow_u[] = area_snow_u[]
 
     δ_d_s = 0
   end
 
   δ_d_s = max(0, δ_d_s)
-  ms_g[] = max(0, ms_g_last + δ_d_s * ρ_new_snow)
+  ms.g = max(0, ms_g_last + δ_d_s * ρ_new_snow)
 
   if δ_d_s > 0
     ρ_snow[] = (ρ_snow[] * d_s + ρ_new_snow * δ_d_s) / (d_s + δ_d_s)
   else
     ρ_snow[] = (ρ_snow[] - 250) * exp(-0.001 * kstep / 3600.0) + 250.0
   end
 
-  d_s = ms_g[] > 0 ? ms_g[] / ρ_snow[] : 0.0
-  perc_snow_g = min(ms_g[] / (0.05 * ρ_snow[]), 1.0)
+  d_s = ms.g > 0 ? ms.g / ρ_snow[] : 0.0
+  perc_snow_g = min(ms.g / (0.05 * ρ_snow[]), 1.0)
 
   if snowrate_o > 0
     albedo_v_snow[] = (albedo_v_snow[] - 0.70) * exp(-0.005 * kstep / 3600) + 0.7
@@ -79,28 +80,30 @@ function snowpack_stage1_jl(temp_air::Float64, prcp::Float64,
     albedo_n_snow[] = albedo_n_Newsnow
   end
 
-  perc_snow_o, perc_snow_u, perc_snow_g, d_s
+  perc_snow = CanopyLayer(perc_snow_o, perc_snow_u, perc_snow_g)
+  perc_snow, d_s
 end
 
 
-function snowpack_stage2_jl(evapo_snow_o::Float64, evapo_snow_u::Float64,
-  ms_o::Ref{Float64}, ms_u::Ref{Float64})
+function snowpack_stage2_jl(evapo_snow_o::Float64, evapo_snow_u::Float64, ms::CanopyLayer{Float64})
+  # ms_o::Ref{Float64}, ms_u::Ref{Float64})
 
   # kstep::Float64 = kstep  # length of step
-  ms_o[] = max(0, ms_o[] - evapo_snow_o * kstep)
-  ms_u[] = max(0, ms_u[] - evapo_snow_u * kstep)
+  ms.o = max(0, ms.o - evapo_snow_o * kstep)
+  ms.u = max(0, ms.u - evapo_snow_u * kstep)
 end
 
 
 function snowpack_stage3_jl(temp_air::Float64, temp_snow::Float64, temp_snow_last::Float64, ρ_snow::Float64,
   d_s, d_w,
   # d_s::Ref{Float64}, d_w::Ref{Float64}, 
-  ms_g::Ref{Float64})
+  ms::CanopyLayer{Float64})
+  # ms_g::Ref{Float64})
 
   # kstep = kstep  # kstep in model
   # it is assumed sublimation happens before the melting and freezing process
   d_s_sup = d_s  # this snow depth has already considered sublimation
-  ms_sup = ms_g[]
+  ms_sup = ms.g
 
   # parameters
   cp_ice = 2228.261
@@ -152,8 +155,8 @@ function snowpack_stage3_jl(temp_air::Float64, temp_snow::Float64, temp_snow_las
   end
 
   # δ in mass of snow on ground
-  ms_g[] = ms_g[] - ms_melt + mw_frozen
-  ms_g[] = max(0, ms_g[])
+  ms.g = ms.g - ms_melt + mw_frozen
+  ms.g = max(0, ms.g)
 
   # δ of depth in snow
   melt_d_s = ms_melt / ρ_snow

Diff for: src/DataType/CanopyLayer.jl

@@ -0,0 +1,34 @@
+import Base.show
+
+abstract type AbstractLayer{FT} end
+
+"""
+x = CanopyLayer{Float64}()
+x = CanopyLayer{Ref{Float64}}()
+"""
+@with_kw_noshow mutable struct CanopyLayer{FT} <: AbstractLayer{FT}
+  o::FT = FT(0.0) # overlayer
+  u::FT = FT(0.0) # underlayer
+  g::FT = FT(0.0) # ground
+end
+
+# 注意，如果是Ref，将共用相同的地址
+CanopyLayer(o::FT) where {FT} = CanopyLayer{FT}(o, o, o)
+
+CanopyLayer(o::FT, u::FT) where {FT} = CanopyLayer{FT}(; o, u)
+
+
+function Base.show(io::IO, x::AbstractLayer{FT}) where {FT}
+  println(io, typeof(x))
+  names = fieldnames(typeof(x))
+  for i in 1:nfields(x)
+    name = names[i]
+    value = getfield(x, name)
+    T = typeof(value)
+    isa(value, Ref) && (value = value[])
+    println(io, "$name: $T $(value)")
+  end
+end
+
+
+export CanopyLayer

Diff for: src/DataType/DataType.jl

@@ -13,6 +13,7 @@ nzero(n) = tuple(zeros(n)...) # n double zero
 
 include("Constant.jl")
 include("Leaf.jl")
+include("CanopyLayer.jl")
 include("Soil.jl")
 
 include("InterTempVars.jl")

Diff for: src/DataType/InterTempVars.jl

@@ -43,32 +43,32 @@ InterTempLeafs(x0) = InterTempLeafs(; x0)
 export SurfaceMass
 
 @with_kw mutable struct SurfaceMass{FT<:AbstractFloat}
-  ρ_snow       ::FT = 0.0
-  prcp_g       ::FT = 0.0
+  ρ_snow::FT = 0.0
+  prcp_g::FT = 0.0
 
-  depth_water  ::FT = 0.0
-  depth_snow   ::FT = 0.0
+  depth_water::FT = 0.0
+  depth_snow::FT = 0.0
 
-  perc_snow_o  ::FT = 0.0 # Xcs_o
-  perc_snow_u  ::FT = 0.0 # Xcs_u
-  perc_snow_g  ::FT = 0.0 # Xcs_g
+  perc_snow_o::FT = 0.0 # Xcs_o
+  perc_snow_u::FT = 0.0 # Xcs_u
+  perc_snow_g::FT = 0.0 # Xcs_g
 
-  perc_water_o ::FT = 0.0 # Xcl_o
-  perc_water_u ::FT = 0.0 # Xcl_u
+  perc_water_o::FT = 0.0 # Xcl_o
+  perc_water_u::FT = 0.0 # Xcl_u
 
-  area_snow_o  ::FT = 0.0 # Ac_snow_o
-  area_snow_u  ::FT = 0.0 # Ac_snow_u
+  area_snow_o::FT = 0.0 # Ac_snow_o
+  area_snow_u::FT = 0.0 # Ac_snow_u
 
   pre_mass_water_o::FT = 0.0 # Wcl_o
   pre_mass_water_u::FT = 0.0 # Wcl_u
-  
+
   mass_water_o::FT = 0.0 # Wcl_o
   mass_water_u::FT = 0.0 # Wcl_u
   mass_water_g::FT = 0.0 # Wcl_g
 
-  mass_snow_o  ::FT = 0.0 # Wcs_o
-  mass_snow_u  ::FT = 0.0 # Wcs_u
-  mass_snow_g  ::FT = 0.0 # Wcs_g
+  mass_snow_o::FT = 0.0 # Wcs_o
+  mass_snow_u::FT = 0.0 # Wcs_u
+  mass_snow_g::FT = 0.0 # Wcs_g
 
   albedo_v_snow::FT = 0.0
   albedo_n_snow::FT = 0.0
@@ -107,13 +107,13 @@ end
   # area
   Ac_snow_o::Vector{FT} = zeros(MAX_Loop)
   Ac_snow_u::Vector{FT} = zeros(MAX_Loop)
-  
+
   rho_snow::Vector{FT} = zeros(MAX_Loop)
   r_rain_g::Vector{FT} = zeros(MAX_Loop)
-  
+
   alpha_v_sw::Vector{FT} = zeros(MAX_Loop)
   alpha_n_sw::Vector{FT} = zeros(MAX_Loop)
-  
+
   Trans_o::Vector{FT} = zeros(MAX_Loop)
   Trans_u::Vector{FT} = zeros(MAX_Loop)
   Eil_o::Vector{FT} = zeros(MAX_Loop)
@@ -144,18 +144,18 @@ function init_vars!(x::InterTempVars)
   # x.Qhc_o .= 0.
   # x.Qhc_u .= 0.
   # x.Qhg .= 0.
-  x.Wcl_o .= 0.
+  x.Wcl_o .= 0.0
   x.Wcl_u .= 0.0
 
   x.Wcs_o .= 0.0
   x.Wcs_u .= 0.0
   x.Wcs_g .= 0.0
 
-  x.Xcl_o .= 0.  
-  x.Xcl_u .= 0.
+  x.Xcl_o .= 0.0
+  x.Xcl_u .= 0.0
 
   x.Xcs_o .= 0.0
-  x.Xcs_u .= 0.
+  x.Xcs_u .= 0.0
   x.Xcs_g .= 0.0
 
   x.rho_snow .= 0.0

Diff for: src/DataType/Leaf.jl

@@ -1,10 +1,4 @@
-export @reset, reset!
-macro reset(x)
-  quote
-    reset!($(esc(x)))
-  end
-end
-
+# they also abstract vector but with name
 @with_kw mutable struct Leaf
   o_sunlit::Cdouble = 0.0
   o_shaded::Cdouble = 0.0
@@ -19,7 +13,6 @@ Leaf(o::Cdouble, u::Cdouble) = Leaf(o, u)
 Base_ops = ((:Base, :+), (:Base, :-), (:Base, :*), (:Base, :/))
 
 for (m, f) in Base_ops
-  # _f = Symbol(m, ".:", f)
   @eval begin
     $m.$f(a::Leaf, b::Leaf) = begin
       Leaf(
@@ -33,31 +26,34 @@ for (m, f) in Base_ops
 end
 
 
+export @reset, reset!
+macro reset(x)
+  quote
+    reset!($(esc(x)))
+  end
+end
+
 reset!(x::Leaf) = init_leaf_dbl(x, 0.0)
 
 function init_leaf_struct(x::Leaf, replacement::Leaf)
   x.o_sunlit = replacement.o_sunlit
   x.o_shaded = replacement.o_shaded
   x.u_sunlit = replacement.u_sunlit
   x.u_shaded = replacement.u_shaded
-  # ccall((:init_leaf_struct, libbeps), Cvoid, (Ptr{Leaf}, Ptr{Leaf}), Ref(x), Ref(replacement))
 end
 
 function init_leaf_dbl(x::Leaf, replacement::Float64)
   x.o_sunlit = replacement
   x.o_shaded = replacement
   x.u_sunlit = replacement
   x.u_shaded = replacement
-  nothing
-  # ccall((:init_leaf_dbl, libbeps), Cvoid, (Ptr{Leaf}, Cdouble), Ref(x), replacement)
 end
 
-function init_leaf_dbl2(x, overstory, understory)
+function init_leaf_dbl2(x::Leaf, overstory, understory)
   x.o_sunlit = overstory
   x.o_shaded = overstory
   x.u_sunlit = understory
   x.u_shaded = understory
-  # ccall((:init_leaf_dbl2, libbeps), Cvoid, (Ptr{Leaf}, Cdouble, Cdouble), Ref(x), overstory, understory)
 end