add PET_1948

example/ex01_栾城.jl

@@ -43,6 +43,7 @@ Tas = cumsum(Tas)
 
 Gi = 0.4 .* Rn .* exp.(-0.5 .* LAI) # G_soil
 s_VODi = (VOD ./ nanmaximum(VOD)) .^ 0.5 # VOD-stress
+topt = 24.0
 
 ET, Tr, Es, Ei, Esb, SM, RF, GW = 
   SiTHv2_site(Rn, Tavg, Tas, Prcp, Pa, Gi, LAI, s_VODi, topt, soilpar, pftpar, state, false)

src/SPAC.jl

@@ -10,6 +10,7 @@ export ZM
 
 # Set the soil depth for three soil layers
 const ZM = [50.0, 1450.0, 3500.0]  # mm
+const atm = 101.325  # atmospheric pressure (kPa)
 
 include("DataType.jl")
 include("Param/get_pftpar.jl")

src/SiTHv2.jl

@@ -31,7 +31,9 @@ export SiTHv2_site
 """
 function SiTHv2(Rn, Ta, Tas, Topt, P, Pa, s_VOD, G, LAI, soilpar, pftpar, state::State; Kc=1.0)
   (; wa, zgw, snowpack) = state
-  pEc, pEs = potentialET(Rn, G, LAI, Ta, Pa; Kc) # PET allocated to canopy and soil surface
+  Kc = [1.0, 1.0, 1.0]
+  VPD, U2, doy = 0.0, 0.0, 0
+  pEc, pEs, ET0 = potentialET(Rn, G, LAI, Ta, Pa, VPD, U2, doy; Kc) # PET allocated to canopy and soil surface
   Ei, fwet, PE = interception(P, pEc, LAI, pftpar)  # Interception evaporation
 
   # Snow sublimation, snow melt
@@ -56,7 +58,7 @@ function _run_model(Rn::T, Ta::T, Tas::T, Prcp::T, Pa::T, G::T, LAI::T, s_VOD::T
   for i in 1:ntime
     _Kc = 180 <= i <= 220 ? Kc : 1.0
     Et, Tr, Es, Ei, Esb, srf, _, _, _ = SiTHv2(Rn[i], Ta[i], Tas[i], Top, Prcp[i], Pa[i], s_VOD[i], G[i], LAI[i], soilpar, pftpar, state; Kc=_Kc)
-    
+
     output[:ETs][i] = Et
     output[:Trs][i] = Tr
     output[:Ess][i] = Es

src/potentialET.jl

@@ -1,52 +1,109 @@
-const Cp = 1013  # Specific heat (J kg-1 C-1)
+# λ: [J/kg] change to [MJ kg-1]
+const Cp = 1.013 * 1e-3  # Specific heat (MJ kg-1 C-1)
 const ϵ = 0.622  # e (unitless) is the ratio of molecular weight of water to dry air (equal to 0.622)
 
+period_wheat(doy::Int) = 32 <= doy <= 166
+period_corn(doy::Int) = 196 <= doy <= 288
+
+# [wheat, corn, non-gw]
+function iGrowthPeriod(doy::Int)
+  i = 3
+  period_wheat(doy) && (i = 1)
+  period_corn(doy) && (i = 2)
+  return i
+end
+
+
 """
 Potential ET partition
 
 ## INPUT:
-- Ta    : air temperature, C
-- Rn    : average daily net radiation, W/m^2
-- Pa    : atmospheric pressure, kPa
-- LAI   : Leaf area index, 1
-- G     : Soil heat flux, W/m^2
+- Ta  : air temperature, C
+- Rn  : average daily net radiation, W/m^2
+- Pa  : atmospheric pressure, kPa
+- LAI : Leaf area index, 1
+- G   : Soil heat flux, W/m^2
+
+## Optional:
+- k   : 消光系数，default `[0.6, 0.6, 0.6]`
+- Kc  : 作物折腾转换系数，default `[0.75, 0.9, 1.0]`
 
 ## OUTPUT
 - pEc   : potential Transpiration, mm/day
 - pEs   : potential Soil evaporation, mm/day
 """
-function potentialET(Rn::T, G::T, LAI::T, Ta::T, Pa::T; Kc::T = 1.0) where {T<:Real}
-  k = 0.6  # the empirical extinction coefficient set as 0.6
+function potentialET(Rn::T, G::T, LAI::T, Ta::T, Pa::T, 
+  VPD::T, U2::T, doy::Int;
+  method="PT1972",
+  Kc::Vector=[0.75, 0.9, 1.0],
+  k::Vector=[0.6, 0.6, 0.6]) where {T<:Real}
+
+  i = iGrowthPeriod(doy)
+  _Kc = Kc[i]
+  _k = k[i]
+
   # Radiation located into soil and canopy, separately
-  Rns = exp(-k * LAI) * Rn
+  Rns = exp(-_k * LAI) * Rn
   Rnc = Rn - Rns
 
-  λ = cal_lambda(Ta)                   # [J/kg]  
-  Δ = cal_slope(Ta)                    # [kPa/degC]
-  γ = (Cp * Pa) / (ϵ * λ)              # Psychrometric constant (kPa/degC)
+  ## Potential Transpiration and Soil evaporation, mm/day
+  if method == "PT1972"
+    ET0 = ET0_PT1972(Rn, Ta, Pa)
+    pEc = ET0_PT1972(Rnc, Ta, Pa) # * _Kc
+    pEs = ET0_PT1972(Rns - G, Ta, Pa)
+  elseif method == "Penman48"
+    ET0 = ET0_Penman48(Rn, Ta, VPD, U2)             # all
+    pEc = ET0_Penman48(Rnc, Ta, VPD, U2, Pa) * _Kc  # canopy
+    pEs = ET0_Penman48(Rns - G, Ta, VPD, U2, Pa)    # soil
+  end
+  return pEc, pEs, ET0
+end
+
 
-  # Potential Transpiration and Soil evaporation, mm/day
-  pEc = ET0_PT1972(Rnc, Δ, γ, λ) * Kc
-  pEs = ET0_PT1972(Rns - G, Δ, γ, λ)
-  return pEc, pEs
+function ET0_eq(Rn::T, Tair::T, Pa::T=atm, args...) where {T<:Real}
+  λ::T = cal_lambda(Tair)     # [MJ kg-1]
+  Δ::T = cal_slope(Tair)      # [kPa degC-1]
+  γ::T = Cp * Pa / (ϵ * λ)    # [kPa degC-1], Psychrometric constant
+
+  Eeq::T = Δ / (Δ + γ) * Rn |> x -> W2mm(x, λ)
+  Eeq = max(Eeq, 0)
+  (; Eeq, λ, Δ, γ)
 end
 
+# - Rn: W m-2
+# - α: PT coefficient for water saturated surface
+function ET0_PT1972(Rn::T, Tair::T, Pa::T=atm; α=1.26) where {T<:Real}
+  Eeq = ET0_eq(Rn, Tair, Pa)[1]
+  α * Eeq # [mm d-1]
+end
 
-# Rn: W m-2
-function ET0_PT1972(Rn::T, Δ::T, γ::T, λ::T) where {T<:Real}
-  α = 1.26  # PT coefficient for water saturated surface
-  ET0 = α * Rn * Δ / (Δ + γ)
-  ET0 = max(ET0, 0)
-  W2mm(ET0, λ)
+"""
+ET0_Penman48(Rn, Tair, VPD, Uz)
+"""
+function ET0_Penman48(Rn::T, Tair::T, VPD::T, Uz::T, Pa::T=atm; z_wind=2) where {T<:Real}
+  Eeq, λ, Δ, γ = ET0_eq(Rn, Tair, Pa)
+  U2::T = cal_U2(Uz, z_wind)
+  Evp::T = γ / (Δ + γ) * 6.43 * (1 + 0.536 * U2) * VPD / λ
+  ET0::T = Eeq + Evp
+  ET0
 end
 
+
 # λ: latent heat of vaporization (J kg-1)
-W2mm(w::T, λ::T) where {T<:Real} = w * 86400 / λ
+W2mm(w::T, λ::T) where {T<:Real} = w * 0.086400 / λ
 
 # λ: latent heat of vaporization (J kg-1)
-cal_lambda(Ta::T) where {T<:Real} = 2.501e6 - 2361 * Ta
+cal_lambda(Ta::T) where {T<:Real} = 2.501 - 2.361 / 1000 * Ta
 
 # cal_gamma(Pa, λ) = (Cp * Pa) / (ϵ * λ)
 cal_es(Ta::T) where {T<:Real} = 0.6108 * exp((17.27 * Ta) / (Ta + 237.3))
 
 cal_slope(Ta::T) where {T<:Real} = 4098 * cal_es(Ta) / ((Ta + 237.3)^2)
+
+function cal_U2(Uz::T, z=10.0) where {T<:Real}
+  z == 2 && (return Uz)
+  Uz * 4.87 / log(67.8 * z - 5.42)
+end
+
+
+export ET0_eq, ET0_Penman48, ET0_PT1972

test/runtests.jl

@@ -1,13 +1,14 @@
 using Test, SPAC
-include("test_栾城_2010.jl")
 
 @testset "potentialET" begin
   Rn = 100.0
   G = 5.0
   LAI = 3.0
   Ta = 20.0
   Pa = 100.0
-  pEc, pEs = potentialET(Rn, G, LAI, Ta, Pa)
+
+  VPD, U2, doy = 0.0, 0.0, 0
+  pEc, pEs, ET0 = potentialET(Rn, G, LAI, Ta, Pa, VPD, U2, doy)
 
   @test pEc ≈ 2.5389604966643824
   @test pEs ≈ 0.3507115521629298
@@ -32,3 +33,5 @@ end
   r = pTr_partition(pEc, fwet, wa1, wa2, wa3, soilpar, pftpar, ZM)
   @test r == (0.004508186497043249, 9.864271368015835, 0.1312204454871218)
 end
+include("test-PET.jl")
+include("test_栾城_2010.jl")

test/test-PET.jl

@@ -0,0 +1,11 @@
+using SPAC, Test
+
+@testset "ET0_Penman48" begin
+  Rn = 200.0    # W/m²
+  Ta = 25.0     # °C
+  VPD = 1.5     # kPa
+  Uz = 3.0      # m/s
+
+  @test ET0_Penman48(Rn, Ta, VPD, Uz) == 7.9265544809634365
+  @test ET0_PT1972(Rn, Ta) == 6.564860193238437
+end