add more radiation function

src/Units/test-Rn.jl

@@ -0,0 +1,10 @@
+## 使用Unit速度会有较大牺牲
+# function cal_Rn(Rs::Quantity{T}, Rln::Quantity{T}, Tavg::Quantity{T},
+#   α::Float64, ϵ::Float64) where {T<:Real}
+#   K = uconvert(u"K", Tavg)
+#   RLout = Stefan * K^4 |> unit(Rs)
+#   Rnl = ϵ * (Rln - RLout)
+#   Rns = (1.0 - α) * Rs
+#   Rn = Rns + Rnl
+#   max(Rn, 0.0 * unit(Rn))
+# end

src/cal_humidity.jl

@@ -51,3 +51,5 @@ function q2VPD(q, Tmin, Tmax, Pa=atm)
   es = cal_es(Tmin, Tmax)
   max(es - ea, 0.0)
 end
+
+export cal_es, Tdew2RH, Tdew2VPD, ea2q, q2ea, q2RH, q2VPD

src/cal_radiation.jl

@@ -1,5 +1,9 @@
+export cal_Rn, cal_Rsi, cal_Rsi_toa, 
+  cal_Rln, cal_Rln_out, cal_Rli, cal_Rln_yang2019
+
+
 """
-    get_Rn(Rs, Rln, Tavg, albedo, emiss)
+    cal_Rn(Rs, Rln, Tavg, albedo, emiss)
 
 Calculate net radiation (Rn) using input parameters.
 
@@ -21,22 +25,194 @@ function cal_Rn(Rs::T, Rln::T, Tavg::T, α::Float64, ϵ::Float64) where {T<:Real
   Rns = (1.0 - α) * Rs
 
   Rn = Rns + Rnl
-  max(Rn, 0.0)
+  Rn
+  # Rn, Rns, Rnl
+  # max(Rn, T(0.0))
+end
+
+# convert from MJ m-2 d-1 to W/m2
+cal_Rln_out(T::FT, ϵ=0.96) where {FT<:Real} = Stefan * (T+K0)^4 / 0.0864 
+
+
+"""
+  cal_Rsi_toa(lat, J)
+
+Estimate daily extraterrestrial radiation in MJ m-2 day-1.
+
+# Arguments
+- `lat`: Latitude in degrees.
+- `J`: Day of the year.
+
+# Returns
+- Extraterrestrial radiation in `MJ m-2 day-1`.
+"""
+function cal_Rsi_toa(lat=0, J::Integer=1)
+  dr = 1 + 0.033 * cos(π * J / 182.5) # Allen, Eq. 23
+  σ = 0.409 * sin(π * J / 182.5 - 1.39) # Allen, Eq. 24
+
+  ws = HourAngleSunSet(lat, J)
+  # 24 * 60 * 0.082 = 118.08
+  lat = deg2rad(lat)
+  Rsi_toa = 118.08 * dr / π * (ws * sin(lat) * sin(σ) + cos(lat) * cos(σ) * sin(ws)) # Allen, Eq. 21
+  max(Rsi_toa, 0)
+end
+
+
+"""
+  cal_Rsi(lat, J, ssd=nothing, cld=nothing, Z=0, a=0.25, b=0.5)
+
+Daily inward shortwave solar radiation at crop surface in MJ m-2 day-1 by
+providing sunshine duration (SSD) in hours or cloud cover in fraction.
+
+# Arguments
+- `lat`: Latitude in degrees.
+- `J`: Day of the year.
+- `ssd`: Sunshine duration in hours. If `ssd` is `nothing`, `Rsi` is the
+  clear-sky solar radiation. Default is `nothing`.
+- `cld`: Cloud cover in fraction. If provided it would be directly used to
+  calculate solar radiation rather than SSD and parameter a and b. Default is
+  `nothing`.
+- `Z`: Elevation in meters. Default is 0.
+- `a`: Coefficient of the Angstrom formula. Determine the relationship between
+  ssd and radiation. Default is 0.25.
+- `b`: Coefficient of the Angstrom formula. Default is 0.50.
+
+# Returns
+- A tuple of solar radiation at crop surface in MJ m-2 day-1:
+  - `Rsi`: Surface downward shortwave radiation.
+  - `Rsi_o`: Clear-sky surface downward shortwave radiation.
+  - `Rsi_toa`: Extraterrestrial radiation.
+"""
+function cal_Rsi(lat=0, J::Integer=1, ssd=nothing, cld=nothing; 
+  Z=0, a=0.25, b=0.5)
+
+  Rsi_toa = cal_Rsi_toa(lat, J)
+
+  if cld !== nothing
+    nN = (1 - cld)
+  else
+    N = SunshineDuration(lat, J)
+    ssd = ssd === nothing ? N : ssd
+    nN = ssd / N
+  end
+
+  Rsi = (a + b * nN) * Rsi_toa
+  Rsi_o = (0.75 + 2 * Z / 1e5) * Rsi_toa
+
+  Rsi, Rsi_o, Rsi_toa
 end
 
-## 使用Unit速度会有较大牺牲
-# function cal_Rn(Rs::Quantity{T}, Rln::Quantity{T}, Tavg::Quantity{T},
-#   α::Float64, ϵ::Float64) where {T<:Real}
 
-#   K = uconvert(u"K", Tavg)
-#   RLout = Stefan * K^4 |> unit(Rs)
+"""
+    cal_Rln(Tmax, Tmin, ea, cld)
+    cal_Rln(Tmax, Tmin, ea, Rsi, Rsi_o)
+
+Net outgoing longwave radiation.
+
+# Arguments
+- `Tmax`: Daily maximum air temperature at 2m height in degrees Celsius.
+- `Tmin`: Daily minimum air temperature at 2m height in degrees Celsius.
+- `ea`: Actual vapor pressure in kPa. Can be estimated by maximum or minimum air
+  temperature and mean relative humidity.
+- `Rsi`: Incoming shortwave radiation at crop surface in MJ m-2 day-1. Default
+  is `nothing`.
+- `Rsi_o`: Clear sky incoming shortwave radiation, i.e., extraterrestrial
+  radiation multiply by clear sky transmissivity (i.e., a + b, a and b are
+  coefficients of Angstrom formula. Normally 0.75) in MJ m-2 day-1. Default is
+  `nothing`.
+- `cld`: Cloud cover in fraction. If provided it would be directly used to
+  calculate net outgoing longwave radiation than Rso. Default is `nothing`.
+
+# Returns
+- Net outgoing longwave radiation in MJ m-2 day-1.
+"""
+function cal_Rln(Tmax, Tmin, ea, cld)
+  cld = isnan(cld) ? 1.0 : cld
+
+  return (4.093e-9 * (((Tmax + 273.15)^4 + (Tmin + 273.15)^4) / 2)) *
+         (0.34 - (0.14 * sqrt(ea))) *
+         (1.35 * (1.0 - cld) - 0.35)
+end
+
+function cal_Rln(Tmax, Tmin, ea, Rsi, Rsi_o)
+  cld = 1.0 - Rsi / Rsi_o
+  cal_Rln(Tmax, Tmin, ea, cld)
+end
+
+
+"""
+    cal_Rln_yang2019(Tₛ, Rsi, Rsi_toa; lat=30, ϵ::Float64=0.96, n1=2.52, n2=2.37, n3=0.035)
 
-#   Rnl = ϵ * (Rln - RLout)
-#   Rns = (1.0 - α) * Rs
-#   Rn = Rns + Rnl
+Calculate the longwave radiation based on the Yang et al. (2019) method.
 
-#   max(Rn, 0.0 * unit(Rn))
-# end
+# Arguments
+- `Ts`: Land surface temperature.
+- `Rsi`: Surface incoming short-wave radiation.
+- `Rsi_toa`: Incoming short-wave radiation at the top of the atmosphere.
+- `lat`: Latitude (in degrees). Default is 30.
+- `ϵ`: Emissivity. Default is 0.96.
+- `n1`, `n2`, `n3`: Parameters for `ΔT`. See Yang 2019 for details.
+
+# References
+1. Yang, Y., & Roderick, M. L. (2019). Radiation, surface temperature and
+   evaporation over wet surfaces. Quarterly Journal of the Royal Meteorological
+   Society, 145(720), 1118–1129. https://doi.org/10.1002/qj.3481
+2. Tu, Z., & Yang, Y. (2022). On the Estimation of Potential Evaporation Under
+   Wet and Dry Conditions. Water Resources Research, 58(4).
+   https://doi.org/10.1029/2021wr031486
+"""
+function cal_Rln_yang2019(Ts, Rsi, Rsi_toa;
+  lat=30, ϵ::Float64=0.96, n1=2.52, n2=2.37, n3=0.035)
+
+  tau = Rsi / Rsi_toa
+  ΔT = n1 * exp(n2 * tau) + n3 * abs(lat)
+  return ϵ * σ * ((Ts - ΔT)^4 - Ts^4)
+end
+
+
+"""
+  cal_Rli(Ta, ea=0, s=1, method="MAR")
+
+Estimate Incoming longwave radiation. Not included in FAO56 paper but added for
+convenience.
 
+# Arguments
+- `Ta`: Near surface air temperature in degrees Celsius.
+- `ea`: Near surface actual vapour pressure in kPa. Default is 0.
+- `s`: Cloud air emissivity, the ratio between actual incoming shortwave
+  radiation and clear sky incoming shortwave radiation. Default is 1.
+- `method`: The method to estimate the air emissivity. Must be one of 'MAR',
+  'SWI', 'IJ', 'BRU', 'SAT', 'KON'. Default is 'MAR'.
 
-export cal_Rn
+# Returns
+- Incoming longwave radiation in W/m².
+
+# References
+1. Satterlund, D. R. (1979), An improved equation for estimating long-wave
+   radiation from the atmosphere, Water Resour. Res., 15( 6), 1649– 1650,
+   doi:10.1029/WR015i006p01649.
+2. Sedlar, J., & Hock, R. (2009). Testing longwave radiation parameterizations
+   under clear and overcast skies at Storglaciären, Sweden. The Cryosphere,
+   3(1), 75-84.
+"""
+function cal_Rli(Ta, ea=0; s=1, method="MAR")
+  ea *= 10 # to hPa
+  Ta += 273.15
+  if method == "MAR"
+    ep_ac = 0.5893 + 5.351e-2 * sqrt(ea / 10)
+  elseif method == "SWI"
+    ep_ac = 9.294e-6 * Ta^2
+  elseif method == "IJ"
+    ep_ac = 1 - 0.261 * exp(-7.77e-4 * (273 - Ta)^2)
+  elseif method == "BRU"
+    ep_ac = 1.24 * (ea / Ta)^(1 / 7)
+  elseif method == "SAT"
+    ep_ac = 1.08 * (1 - exp(-ea^(Ta / 2016)))
+  elseif method == "KON"
+    a = 0.4393
+    b = 7
+    ep_ac = 0.23 + a * (ea / 10 / Ta)^(1 / b)
+  end
+  ep_a = 1 - s + s * ep_ac
+  return ep_a * σ * Ta^4
+end

src/cal_sun_angle.jl

@@ -1,7 +1,3 @@
-deg2rad(x) = x / 180 * pi
-rad2deg(x) = x / pi * 180
-
-
 function get_md(J)
   date_origin = DateTime(2010, 1, 1, 0, 0, 0) + Day(J - 1)
   Dates.format(date_origin, "mm-dd")
@@ -20,30 +16,43 @@ end
 
 
 """
-    get_σ(J; to_deg=false)
+    SolarDeclinationAngle(J; to_deg=false)
+
+# Arguments
+- σ: Solar Declination Angle, 黄赤交角（太阳赤纬角）
+- ϕ: `纬度`
+- ω: `时角`
+"""
+function SolarDeclinationAngle(J::Integer; deg=false)
+  σ = 0.409 * sin(2 * pi / 365 * J .- 1.39) # Allen, Eq. 24
+  deg ? rad2deg(σ) : σ
+end
 
+"""
+  HourAngleSunSet(lat, J)
 
-- `黄赤交角` : σ, Solar Declination Angle
+Calculate the sunset hour angle.
 
-- `纬度`: ϕ
+# Arguments
+- `lat`: Latitude in degrees.
+- `J`: Day of the year.
 
-- `时角`: ω
+# Returns
+- Sunset hour angle in radian
 """
-function get_σ(J; to_deg=false)
-  σ = 0.409 .* sin.(2 * pi / 365 * J .- 1.39) # in [rad]
-  if to_deg
-    σ = rad2deg.(σ)
-  end
-  σ
+function HourAngleSunSet(lat=0, J::Integer=1; deg=false)
+  lat = deg2rad(lat)
+  σ = SolarDeclinationAngle(J)
+
+  tmp = clamp(-tan(lat) * tan(σ), -1, 1)
+  ws = acos(tmp) # Eq. 25
+  deg ? rad2deg(ws) : ws
 end
 
-function get_ssh(ϕ, J)
-  σ = get_σ.(J)
-  # rad2deg(σ)
-  tmp = -tan.(ϕ) .* tan.(σ)
-  clamp!(tmp, -1, 1)
-  # constrain in the range of [-1, 1]
-  w = acos.(tmp)
+function SunshineDuration(lat=0, J::Integer=1)
+  w = HourAngleSunSet(lat, J)
   w_hour = w / pi * 12 # 距离中午的时间
   w_hour * 2 # ssh
 end
+
+export HourAngleSunSet, SunshineDuration, ssh2time

src/constant.jl

@@ -46,5 +46,6 @@ Cp = 1.013 * 1e-3 # MJ kg-1 degC-1
 
 Stefan = 4.903e-9 # u"MJ / (K^4 * m^2 *d)" # Stefan-Boltzmann constant
 
+σ = 5.67e-8 # 5.67*1e-8 Wm−2 K−4
 
 export atm, Stefan

src/unit_convert.jl

@@ -62,4 +62,5 @@ deg2rad(deg::Real) = deg / 180.0 * π
 rad2deg(rad::Real) = rad / π * 180.0
 
 
-export R2Q, Q2R, MJ2W, MJ2mm, W2MJ, W2mm, F2C, C2F, K2C, C2K, deg2rad, rad2deg
+export R2Q, Q2R, MJ2W, MJ2mm, W2MJ, W2mm, F2C, C2F, K2C, C2K
+# deg2rad, rad2deg

test/runtests.jl

@@ -48,6 +48,6 @@ end
         (NSE=0.8787878787878788, R2=1.0, KGE=0.8181818181818181, R=1.0, RMSE=1.0, MAE=1.0, bias=1.0, bias_perc=18.181818181818183, n_valid=10)
 end
 
-
+include("test-radiation.jl")
 include("test-PMLV2.jl")
 include("test-sceua.jl")

test/test-radiation.jl

@@ -0,0 +1,19 @@
+using Test
+
+@testset "radiation" begin
+  cal_Rsi()
+  cal_Rsi_toa()
+
+  cal_Rln(30, 20, 1, 0.5)
+  cal_Rln_out(10)
+  cal_Rli(10)
+  cal_Rln_yang2019(20, 100, 200)
+
+
+  Rs = 200
+  Rln = 400
+  Tavg = 20
+  albedo = 0.2
+  emiss = 0.96
+  @test_nowarn cal_Rn(Rs, Rln, Tavg, albedo, emiss)
+end