mv Base to BEPS

Author: kongdd

src/BEPS.jl

@@ -21,6 +21,7 @@ libbeps = path_proj("deps/libbeps.dll")
 # include("DataFrames.jl")
 # include("Ipaper.jl")
 # include("c2julia.jl")
+include("SPAC/SPAC.jl")
 include("DataType/DataType.jl")
 include("BEPS_modules.jl")
 

src/BEPS_modules.jl

@@ -1,4 +1,4 @@
-include("Base/Base.jl")
+# include("Base/Base.jl")
 include("Param/Parameter.jl")
 
 include("Soil/Soil.jl")

src/Base/Base.jl

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

src/Base/helpers.jl

@@ -0,0 +1,16 @@
+
+## BEPS modules
+function update_Gw!(Gw::Leaf, Gs_new::Leaf, Ga_o, Ga_u, Gb_o, Gb_u)
+  Gw.o_sunlit = 1.0 / (1.0 / Ga_o + 1.0 / Gb_o + 1.0 / Gs_new.o_sunlit)  #conductance for water
+  Gw.o_shaded = 1.0 / (1.0 / Ga_o + 1.0 / Gb_o + 1.0 / Gs_new.o_shaded)
+  Gw.u_sunlit = 1.0 / (1.0 / Ga_u + 1.0 / Gb_u + 1.0 / Gs_new.u_sunlit)
+  Gw.u_shaded = 1.0 / (1.0 / Ga_u + 1.0 / Gb_u + 1.0 / Gs_new.u_shaded)
+end
+
+function update_Gc!(Gc::Leaf, Gs_new::Leaf, Ga_o, Ga_u, Gb_o, Gb_u)
+  Gc.o_sunlit = 1.0 / (1.0 / Ga_o + 1.4 / Gb_o + 1.6 / Gs_new.o_sunlit) # conductance for CO2
+  Gc.o_shaded = 1.0 / (1.0 / Ga_o + 1.4 / Gb_o + 1.6 / Gs_new.o_shaded)
+  Gc.u_sunlit = 1.0 / (1.0 / Ga_u + 1.4 / Gb_u + 1.6 / Gs_new.u_sunlit)
+  Gc.u_shaded = 1.0 / (1.0 / Ga_u + 1.4 / Gb_u + 1.6 / Gs_new.u_shaded)
+end
+

src/DataType/DataType.jl

@@ -63,18 +63,3 @@ function multiply!(Z::Leaf, X::Leaf, Y::Leaf)
   Z.u_sunlit = X.u_sunlit * Y.u_sunlit
   Z.u_shaded = X.u_shaded * Y.u_shaded
 end
-
-
-function update_Gw!(Gw::Leaf, Gs_new::Leaf, Ga_o, Ga_u, Gb_o, Gb_u)
-  Gw.o_sunlit = 1.0 / (1.0 / Ga_o + 1.0 / Gb_o + 1.0 / Gs_new.o_sunlit)  #conductance for water
-  Gw.o_shaded = 1.0 / (1.0 / Ga_o + 1.0 / Gb_o + 1.0 / Gs_new.o_shaded)
-  Gw.u_sunlit = 1.0 / (1.0 / Ga_u + 1.0 / Gb_u + 1.0 / Gs_new.u_sunlit)
-  Gw.u_shaded = 1.0 / (1.0 / Ga_u + 1.0 / Gb_u + 1.0 / Gs_new.u_shaded)
-end
-
-function update_Gc!(Gc::Leaf, Gs_new::Leaf, Ga_o, Ga_u, Gb_o, Gb_u)
-  Gc.o_sunlit = 1.0 / (1.0 / Ga_o + 1.4 / Gb_o + 1.6 / Gs_new.o_sunlit) # conductance for CO2
-  Gc.o_shaded = 1.0 / (1.0 / Ga_o + 1.4 / Gb_o + 1.6 / Gs_new.o_shaded)
-  Gc.u_sunlit = 1.0 / (1.0 / Ga_u + 1.4 / Gb_u + 1.6 / Gs_new.u_sunlit)
-  Gc.u_shaded = 1.0 / (1.0 / Ga_u + 1.4 / Gb_u + 1.6 / Gs_new.u_shaded)
-end

src/SPAC/BEPS_helper.jl

@@ -0,0 +1,60 @@
+using DocStringExtensions: TYPEDFIELDS
+using DataFrames: DataFrame
+
+# include("../SPAC/SPAC.jl")
+
+export s_coszs, lai2, VCmax
+export fill_meteo!, meteo_pack_jl
+export snow_density
+
+include("Leaf.jl")
+include("list.jl")
+include("lai2.jl")
+include("VCmax.jl")
+include("snow_density.jl")
+include("helper.jl")
+include("BEPS_helper.jl")
+
+"""
+    s_coszs(jday::Int, j::Int, lat::Float64, lon::Float64)
+
+# Example
+```julia
+jday, hour, lat, lon = 20, 12, 20., 120.
+s_coszs(jday, hour, lat, lon)
+```
+"""
+function s_coszs(jday::Int, hour::Int, lat::Float64, lon::Float64)
+  Delta = 0.006918 - 0.399912 * cos(jday * 2π / 365.0) + 0.070257 * sin(jday * 2π / 365.0) -
+          0.006758 * cos(jday * 4π / 365.0) + 0.000907 * sin(jday * 4π / 365.0)
+  # delta is the declination angle of sun.
+
+  hr = hour + lon / 15.0  # UTC time
+  # hr =j*24.0/RTIMES; # local time
+  hr > 24 && (hr = hr - 24)
+  hr < 0 && (hr = 24 + hr)
+
+  Lat_arc = π * lat / 180.0
+  Hsolar1 = (hr - 12.0) * 2.0 * π / 24.0 # local hour angle in arc.
+
+  # sin(h)
+  CosZs = cos(Delta) * cos(Lat_arc) * cos(Hsolar1) + sin(Delta) * sin(Lat_arc)
+  return CosZs
+end
+
+function meteo_pack_jl(Ta::FT, RH::FT) where {FT<:Real}
+  ρₐ::FT = 1.292 # ρ_air, kg/m3
+  es::FT = cal_es(Ta)
+  ea::FT = es * RH / 100
+  VPD::FT = es - ea
+
+  q::FT = ea2q(ea)
+  cp::FT = cal_cp(q)
+
+  λ::FT = cal_lambda(Ta)
+  Δ::FT = cal_slope(Ta) # slope of es
+  γ::FT = 0.066         # kPa/K, 
+  # lambda = cal_lambda(Ta) # J kg-1
+  # psy = cp * 101.13 / (0.622 * lambda)  
+  (; ρₐ, cp, VPD, λ, Δ, γ, es, ea, q)
+end

src/DataType/Leaf.jl renamed to src/SPAC/Leaf.jl

@@ -3,7 +3,7 @@
 # VCmax25 = 62.5          # maximum capacity of Rubisco at 25C-VCmax	
 # N_leaf = 3.10 + 1.35    # leaf Nitrogen content	mean value + 1 SD g/m2 
 # slope = 20.72 / 62.5       # slope of VCmax-N curve
-function VCmax(lai::FT, Ω::FT, CosZs::FT, VCmax25::FT, N_leaf::FT, slope::FT)
+function VCmax(lai::FT, Ω::FT, CosZs::FT, VCmax25::FT, N_leaf::FT, slope::FT) where {FT<:Real}
   CosZs <= 0 && return 0.0, 0.0 # 光合仅发生在白天
   K = 0.5 * Ω / CosZs # assuming a spherical leaf angle distribution
   Kn = 0.3            # 0.713/2.4

src/SPAC/SPAC.jl

@@ -0,0 +1,67 @@
+# pow = ^
+# pow(x::FT, y)::FT where {FT<:Real} = x^y
+pow(x, y) = x^y
+
+function blackbody(T::FT) where {FT<:Real}
+  σ = 5.67 / 100000000    # stephen-boltzman constant
+  σ * (T + 273.15)^4
+end
+
+function cal_Rln(emiss::FT, T::FT) where {FT<:Real}
+  emiss * blackbody(T)
+end
+
+# kPa deg-1
+function cal_slope(Ta::FT)::FT where {FT<:Real}
+  2503.0 / pow((Ta + 237.3), 2) * exp(17.27 * Ta / (Ta + 237.3))
+end
+
+# kPa
+function cal_es(Ta::FT)::FT where {FT<:Real} 
+  0.61078 * exp(17.3 * Ta / (237.3 + Ta))
+end
+
+function cal_ea(Ta::FT, RH::FT)::FT where {FT<:Real}
+  cal_es(Ta) * RH / 100
+end
+
+function cal_lambda(Ta::FT)::FT where {FT<:Real} 
+  (2.501 - 0.00237 * Ta) * 1000000
+end
+
+function ea2q(ea::FT, Pa::FT=101.35)::FT where {FT<:Real} 
+  0.622 * ea / (Pa - 0.378 * ea)
+end
+
+function RH2q(Ta::FT, RH::FT)::FT where {FT<:Real}
+  es = cal_es(Ta)
+  ea = es * RH / 100
+  ea2q(ea)
+end
+
+"""
+# Arguments
+- `q`  : specific humidity, g / kg
+- `tem`: air temperature, ℃
+"""
+function q2RH(q::FT, tem::FT)::FT where {FT<:Real}
+  # Vapour pressure in mbar
+  ea = 0.46 * q * (tem + 273.16) / 100
+  es = 6.1078 * exp((17.269 * tem) / (237.3 + tem))
+  clamp(ea / es * 100, 0.0, 100.0)
+end
+
+
+function cal_cp(q::FT)::FT where {FT<:Real} 
+  1004.65 * (1 + 0.84 * q)
+end
+
+function cal_cp(Ta::FT, RH::FT)::FT where {FT<:Real}
+  q = RH2q(Ta, RH)
+  cal_cp(q)
+end
+
+export pow, cal_ea, cal_es,
+  cal_slope, cal_lambda, cal_cp, ea2q, RH2q,
+  cal_Rln,
+  blackbody

src/Base/VCmax.jl renamed to src/SPAC/VCmax.jl

@@ -34,30 +34,3 @@ function partition_lai(lai, Ω, CosZs)
   sunlit, shaded
 end
 
-
-"""
-    s_coszs(jday::Int, j::Int, lat::Float64, lon::Float64)
-
-# Example
-```julia
-jday, hour, lat, lon = 20, 12, 20., 120.
-s_coszs(jday, hour, lat, lon)
-```
-"""
-function s_coszs(jday::Int, hour::Int, lat::Float64, lon::Float64)
-  Delta = 0.006918 - 0.399912 * cos(jday * 2π / 365.0) + 0.070257 * sin(jday * 2π / 365.0) - 
-    0.006758 * cos(jday * 4π / 365.0) + 0.000907 * sin(jday * 4π / 365.0)
-  # delta is the declination angle of sun.
-
-  hr = hour  + lon / 15.0  # UTC time
-  # hr =j*24.0/RTIMES; # local time
-  hr > 24 && (hr = hr - 24)
-  hr < 0 && (hr = 24 + hr)
-
-  Lat_arc = π * lat / 180.0
-  Hsolar1 = (hr - 12.0) * 2.0 * π / 24.0 # local hour angle in arc.
-
-  # sin(h)
-  CosZs = cos(Delta) * cos(Lat_arc) * cos(Hsolar1) + sin(Delta) * sin(Lat_arc)
-  return CosZs
-end