From db6c78c8ce968ff0d425c891a8817d90f7154695 Mon Sep 17 00:00:00 2001 From: Dongdong Kong Date: Tue, 16 Jan 2024 21:15:36 +0800 Subject: [PATCH] separate `cal_Rln_Longwave` --- {scripts => data}/s1_GEE.qmd | 0 src/BEPS/Leaf_Temperature.jl | 20 ++--- src/BEPS/netRadiation.jl | 80 +++++++++++--------- src/BEPS/photosynthesis.jl | 2 +- src/BEPS/surface_temperature.jl | 77 +++++++++---------- src/DataType/InterTempVars.jl | 5 +- src/beps_inter_prg.jl | 127 ++++++++++++-------------------- 7 files changed, 148 insertions(+), 163 deletions(-) rename {scripts => data}/s1_GEE.qmd (100%) diff --git a/scripts/s1_GEE.qmd b/data/s1_GEE.qmd similarity index 100% rename from scripts/s1_GEE.qmd rename to data/s1_GEE.qmd diff --git a/src/BEPS/Leaf_Temperature.jl b/src/BEPS/Leaf_Temperature.jl index bc1cd69..c6ad815 100644 --- a/src/BEPS/Leaf_Temperature.jl +++ b/src/BEPS/Leaf_Temperature.jl @@ -12,18 +12,20 @@ end function Leaf_Temperatures_jl(Tair::Float64, slope::Float64, psychrometer::Float64, VPD_air::Float64, Cp_ca::Float64, Gw::Leaf, Gww::Leaf, Gh::Leaf, - Xcs_o::Float64, Xcl_o::Float64, Xcs_u::Float64, Xcl_u::Float64, + Xcs_o::Float64, Xcl_o::Float64, + Xcs_u::Float64, Xcl_u::Float64, radiation::Leaf, Tc::Leaf) - Tc.o_sunlit = Leaf_Temperature_jl(Tair, slope, psychrometer, VPD_air, Cp_ca, - Gw.o_sunlit, Gww.o_sunlit, Gh.o_sunlit, Xcs_o + Xcl_o, radiation.o_sunlit, true) + args = (Tair, slope, psychrometer, VPD_air, Cp_ca) + Tc.o_sunlit = Leaf_Temperature_jl(args..., + Gw.o_sunlit, Gww.o_sunlit, Gh.o_sunlit, Xcs_o + Xcl_o, radiation.o_sunlit) - Tc.o_shaded = Leaf_Temperature_jl(Tair, slope, psychrometer, VPD_air, Cp_ca, - Gw.o_shaded, Gww.o_shaded, Gh.o_shaded, Xcs_o + Xcl_o, radiation.o_shaded, true) + Tc.o_shaded = Leaf_Temperature_jl(args..., + Gw.o_shaded, Gww.o_shaded, Gh.o_shaded, Xcs_o + Xcl_o, radiation.o_shaded) - Tc.u_sunlit = Leaf_Temperature_jl(Tair, slope, psychrometer, VPD_air, Cp_ca, - Gw.u_sunlit, Gww.u_sunlit, Gh.u_sunlit, Xcs_u + Xcl_u, radiation.u_sunlit, true) + Tc.u_sunlit = Leaf_Temperature_jl(args..., + Gw.u_sunlit, Gww.u_sunlit, Gh.u_sunlit, Xcs_u + Xcl_u, radiation.u_sunlit) - Tc.u_shaded = Leaf_Temperature_jl(Tair, slope, psychrometer, VPD_air, Cp_ca, - Gw.u_shaded, Gww.u_shaded, Gh.u_shaded, Xcs_u + Xcl_u, radiation.u_shaded, true) + Tc.u_shaded = Leaf_Temperature_jl(args..., + Gw.u_shaded, Gww.u_shaded, Gh.u_shaded, Xcs_u + Xcl_u, radiation.u_shaded) end diff --git a/src/BEPS/netRadiation.jl b/src/BEPS/netRadiation.jl index 2fb7bcd..e3dd1fb 100644 --- a/src/BEPS/netRadiation.jl +++ b/src/BEPS/netRadiation.jl @@ -1,4 +1,3 @@ -# @timeit_all function netRadiation_jl(Rs_global::FT, CosZs::FT, temp_o::FT, temp_u::FT, temp_g::FT, lai_o::FT, lai_u::FT, lai_os::FT, lai_us::FT, @@ -29,12 +28,9 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT, α_g::FT = 0.5 * (α_v_gs + α_n_gs) # separate global solar radiation into direct and diffuse one - if (CosZs < 0.001) # solar zenith angle small, all diffuse radiation - ratio_cloud = 0.0 - else - ratio_cloud = Rs_global / (1367 * CosZs) # Luo2018, A4 - end - + # solar zenith angle small, all diffuse radiation + ratio_cloud = (CosZs < 0.001) ? 0.0 : Rs_global / (1367 * CosZs) # Luo2018, A4 + if (ratio_cloud > 0.8) Ra.Rs_df = 0.13 * Rs_global # Luo2018, A2 else @@ -46,7 +42,6 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT, Ra.Rs_dir = Rs_global - Ra.Rs_df # Luo2018, A3 # fraction at each layer of canopy, direct and diffuse. use Leaf only lai here - τ_o_dir::FT = exp(-0.5 * Ω * lai_o / CosZs) τ_u_dir::FT = exp(-0.5 * Ω * lai_u / CosZs) @@ -60,17 +55,8 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT, τ_u_df::FT = exp(-0.5 * Ω * lai_u / cosQ_u) τ_us_df::FT = exp(-0.5 * Ω * lai_us / cosQ_u) - # ϵivity of each part - ea = cal_ea(temp_air, rh) - ϵ_air = 1.0 - exp(-(pow(ea * 10.0, (temp_air + 273.15) / 1200.0))) - ϵ_air = clamp(ϵ_air, 0.7, 1.0) - - ϵ_o = 0.98 - ϵ_u = 0.98 - ϵ_g = 0.96 - + # net short direct radiation on canopy and ground if Rs_global > zero && CosZs > zero - # net short direct radiation on canopy and ground Ra.Rns_o_dir = Ra.Rs_dir * ((1.0 - α_o) - (1.0 - α_u) * τ_o_dir) # dir into dif_under Ra.Rns_u_dir = Ra.Rs_dir * τ_o_dir * ((1.0 - α_u) - (1.0 - α_g) * τ_u_dir) Ra.Rns_g_dir = Ra.Rs_dir * τ_o_dir * τ_u_dir * (1.0 - α_g) @@ -80,8 +66,8 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT, Ra.Rns_g_dir = 0.0 end + # net short diffuse radiation on canopy and ground if Rs_global > zero && CosZs > zero - # net short diffuse radiation on canopy and ground Ra.Rns_o_df = Ra.Rs_df * ((1.0 - α_o) - (1.0 - α_u) * τ_o_df) + 0.21 * Ω * Ra.Rs_dir * (1.1 - 0.1 * lai_o) * exp(-CosZs) # A8 Ra.Rns_u_df = Ra.Rs_df * τ_o_df * ((1.0 - α_u) - (1.0 - α_g) * τ_u_df) + @@ -98,22 +84,7 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT, Rns_u = Ra.Rns_u_dir + Ra.Rns_u_df Rns_g = Ra.Rns_g_dir + Ra.Rns_g_df - # 计算植被和地面的净长波辐射 - Rl_air = cal_Rln(ϵ_air, temp_air) - Rl_o = cal_Rln(ϵ_o, temp_o) - Rl_u = cal_Rln(ϵ_u, temp_u) - Rl_g = cal_Rln(ϵ_g, temp_g) - - Rnl_o = (ϵ_o * (Rl_air + Rl_u * (1.0 - τ_u_df) + Rl_g * τ_u_df) - 2 * Rl_o) * - (1.0 - τ_o_df) + - ϵ_o * (1.0 - ϵ_u) * (1.0 - τ_u_df) * (Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df)) - - Rnl_u = (ϵ_u * (Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df) + Rl_g) - 2 * Rl_u) * (1.0 - τ_u_df) + - (1.0 - ϵ_g) * ((Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df)) * τ_u_df + Rl_u * (1.0 - τ_u_df)) + - ϵ_u * (1.0 - ϵ_o) * (Rl_u * (1.0 - τ_u_df) + Rl_g * τ_u_df) * (1.0 - τ_o_df) - - Rnl_g = ϵ_g * ((Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df)) * τ_u_df + Rl_u * (1.0 - τ_u_df)) - - Rl_g + (1.0 - ϵ_u) * Rl_g * (1.0 - τ_u_df) + Rnl_o, Rnl_u, Rnl_g = cal_Rln_Longwave(temp_air, rh, temp_o, temp_u, temp_g, lai_o, lai_u, Ω) # 计算植被和地面的总净辐射 Rn_o = Rns_o + Rnl_o @@ -161,4 +132,41 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT, end -# TODO: 切分成两个函数，长波与短波 +function cal_Rln_Longwave(temp_air::FT, rh::FT, temp_o::FT, + temp_u::FT, temp_g::FT, lai_o::FT, lai_u::FT, Ω::FT) where {FT<:Real} + + # indicators to describe leaf distribution angles in canopy. slightly related with LAI + cosQ_o::FT = 0.537 + 0.025 * lai_o # Luo2018, A10, a representative zenith angle for diffuse radiation transmission + cosQ_u::FT = 0.537 + 0.025 * lai_u + + τ_o_df::FT = exp(-0.5 * Ω * lai_o / cosQ_o) + τ_u_df::FT = exp(-0.5 * Ω * lai_u / cosQ_u) + + # ϵ of each part + ea = cal_ea(temp_air, rh) + ϵ_air = 1.0 - exp(-(pow(ea * 10.0, (temp_air + 273.15) / 1200.0))) + ϵ_air = clamp(ϵ_air, 0.7, 1.0) + + ϵ_o = 0.98 + ϵ_u = 0.98 + ϵ_g = 0.96 + + # 计算植被和地面的净长波辐射 + Rl_air = cal_Rln(ϵ_air, temp_air) + Rl_o = cal_Rln(ϵ_o, temp_o) + Rl_u = cal_Rln(ϵ_u, temp_u) + Rl_g = cal_Rln(ϵ_g, temp_g) + + Rnl_o = (ϵ_o * (Rl_air + Rl_u * (1.0 - τ_u_df) + Rl_g * τ_u_df) - 2 * Rl_o) * + (1.0 - τ_o_df) + + ϵ_o * (1.0 - ϵ_u) * (1.0 - τ_u_df) * (Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df)) + + Rnl_u = (ϵ_u * (Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df) + Rl_g) - 2 * Rl_u) * (1.0 - τ_u_df) + + (1.0 - ϵ_g) * ((Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df)) * τ_u_df + Rl_u * (1.0 - τ_u_df)) + + ϵ_u * (1.0 - ϵ_o) * (Rl_u * (1.0 - τ_u_df) + Rl_g * τ_u_df) * (1.0 - τ_o_df) + + Rnl_g = ϵ_g * ((Rl_air * τ_o_df + Rl_o * (1.0 - τ_o_df)) * τ_u_df + Rl_u * (1.0 - τ_u_df)) - + Rl_g + (1.0 - ϵ_u) * Rl_g * (1.0 - τ_u_df) + + Rnl_o, Rnl_u, Rnl_g +end diff --git a/src/BEPS/photosynthesis.jl b/src/BEPS/photosynthesis.jl index 32f6a61..e53faec 100644 --- a/src/BEPS/photosynthesis.jl +++ b/src/BEPS/photosynthesis.jl @@ -307,7 +307,7 @@ end rate * exp(tprime * eact / (tref * rugc * t_lk)) end -# + function LAMBDA(tak::Float64)::Float64 y = 3149000.0 - 2370.0 * tak # add heat of fusion for melting ice diff --git a/src/BEPS/surface_temperature.jl b/src/BEPS/surface_temperature.jl index 9094eae..1a3d89c 100644 --- a/src/BEPS/surface_temperature.jl +++ b/src/BEPS/surface_temperature.jl @@ -1,8 +1,9 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_water::FT, capacity_heat_soil1::FT, capacity_heat_soil0::FT, Gheat_g::FT, - depth_soil1::FT, density_snow::FT, tempL_u::FT, netRad_g::FT, - evapo_soil::FT, evapo_water_g::FT, evapo_snow_g::FT, lambda_soil1::FT, - percent_snow_g::FT, heat_flux_soil1::FT, T_ground_last::FT, + depth_soil1::FT, ρ_snow::FT, tempL_u::FT, Rn_g::FT, + evapo_soil::FT, evapo_water_g::FT, evapo_snow_g::FT, λ_soil1::FT, + perc_snow_g::FT, G_soil1::FT, + T_ground_last::FT, T_soil1_last::FT, T_any0_last::FT, T_soil0_last::FT, T_snow_last::FT, T_snow1_last::FT, T_snow2_last::FT) @@ -11,16 +12,16 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat latent_water::FT = cal_lambda(T_air) # J kg-1 latent_snow::FT = 2.83 * 1000000 - density_air::FT = rho_a - cp_air::FT = cal_cp(T_air, rh_air) + ρ_air::FT = rho_a + cp::FT = cal_cp(T_air, rh_air) ra_g::FT = 1.0 / Gheat_g # aerodynamic resistance of heat # thermal conductivity of snow - lambda_snow::FT = 0.021 + 4.2 * density_snow / 10000 + 2.2 * density_snow^3 * 1e-9 + λ_snow::FT = 0.021 + 4.2 * ρ_snow / 10000 + 2.2 * ρ_snow^3 * 1e-9 # available energy on ground for - Gg::FT = netRad_g - evapo_snow_g * latent_snow - (evapo_water_g + evapo_soil) * latent_water + Gg::FT = Rn_g - evapo_snow_g * latent_snow - (evapo_water_g + evapo_soil) * latent_water T_ground::FT = 0.0 T_any0::FT = 0.0 @@ -28,12 +29,12 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat T_snow::FT = 0.0 T_snow1::FT = 0.0 T_snow2::FT = 0.0 - heat_flux::FT = 0.0 + G::FT = 0.0 if depth_snow <= 0.02 ttt = capacity_heat_soil1 * 0.02 / length_step - T_ground = (T_ground_last * ttt * ra_g * depth_soil1 + Gg * ra_g * depth_soil1 + density_air * cp_air * T_air * depth_soil1 + ra_g * lambda_soil1 * T_soil1_last) / - (density_air * cp_air * depth_soil1 + ra_g * lambda_soil1 + ttt * ra_g * depth_soil1) + T_ground = (T_ground_last * ttt * ra_g * depth_soil1 + Gg * ra_g * depth_soil1 + ρ_air * cp * T_air * depth_soil1 + ra_g * λ_soil1 * T_soil1_last) / + (ρ_air * cp * depth_soil1 + ra_g * λ_soil1 + ttt * ra_g * depth_soil1) T_ground = clamp(T_ground, T_ground_last - 25, T_ground_last + 25) T_any0 = T_ground @@ -42,32 +43,32 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat T_snow1 = T_any0 T_snow2 = T_any0 - heat_flux = 2 * lambda_soil1 * (T_any0 - T_soil1_last) / depth_soil1 - heat_flux = clamp(heat_flux, -100.0, 100.0) + G = 2 * λ_soil1 * (T_any0 - T_soil1_last) / depth_soil1 + G = clamp(G, -100.0, 100.0) elseif depth_snow > 0.02 && depth_snow <= 0.05 ttt = capacity_heat_soil1 * 0.02 / length_step # for soil fraction part - T_soil0 = (T_soil0_last * ttt * ra_g * depth_soil1 + Gg * ra_g * depth_soil1 + density_air * cp_air * T_air * depth_soil1 + 2 * ra_g * lambda_soil1 * T_soil1_last) / - (density_air * cp_air * depth_soil1 + 2 * ra_g * lambda_soil1 + ttt * ra_g * depth_soil1) + T_soil0 = (T_soil0_last * ttt * ra_g * depth_soil1 + Gg * ra_g * depth_soil1 + ρ_air * cp * T_air * depth_soil1 + 2 * ra_g * λ_soil1 * T_soil1_last) / + (ρ_air * cp * depth_soil1 + 2 * ra_g * λ_soil1 + ttt * ra_g * depth_soil1) T_soil0 = clamp(T_soil0, T_air - 25.0, T_air + 25.0) - ttt = cp_ice * density_snow * depth_snow / length_step # for snow part - T_snow = (T_snow_last * ttt * ra_g * depth_snow + Gg * ra_g * depth_snow + density_air * cp_air * tempL_u * depth_snow + ra_g * lambda_snow * T_any0_last) / - (density_air * cp_air * depth_snow + ra_g * lambda_snow + ttt * ra_g * depth_snow) + 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 + ρ_air * cp * tempL_u * depth_snow + ra_g * λ_snow * T_any0_last) / + (ρ_air * cp * depth_snow + ra_g * λ_snow + ttt * ra_g * depth_snow) T_snow = clamp(T_snow, T_air - 25.0, T_air + 25.0) - ttt = (lambda_soil1 * T_soil1_last / depth_soil1 + T_snow * lambda_snow + 0.02 * capacity_heat_soil1 / length_step * T_any0_last) / - (lambda_soil1 / depth_soil1 + lambda_snow / depth_snow + 0.02 * capacity_heat_soil1 / length_step) - T_any0 = T_soil0 * (1 - percent_snow_g) + ttt * percent_snow_g + ttt = (λ_soil1 * T_soil1_last / depth_soil1 + T_snow * λ_snow + 0.02 * capacity_heat_soil1 / length_step * T_any0_last) / + (λ_soil1 / depth_soil1 + λ_snow / depth_snow + 0.02 * capacity_heat_soil1 / length_step) + T_any0 = T_soil0 * (1 - perc_snow_g) + ttt * perc_snow_g - heat_flux_snow = lambda_snow / (depth_snow + 0.5 * depth_soil1) * (T_snow - T_soil1_last) - heat_flux_soil = heat_flux_snow * (T_any0 - T_soil1_last) / depth_soil1 + 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 - heat_flux = heat_flux_snow * percent_snow_g + heat_flux_soil * (1 - percent_snow_g) - heat_flux = clamp(heat_flux, -100.0, 100.0) + 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 > zero && T_snow_last <= zero && depth_snow > zero @@ -79,30 +80,30 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat T_snow = 0.0 end - # percent_snow_g =min(1.0,Wg_snow[kkk] / (0.05 * rho_snow[kkk])); # use the fraction before - T_ground = T_snow * percent_snow_g + T_soil0 * (1 - percent_snow_g) + # perc_snow_g =min(1.0,Wg_snow[kkk] / (0.05 * rho_snow[kkk])); # use the fraction before + T_ground = T_snow * perc_snow_g + T_soil0 * (1 - perc_snow_g) T_ground = clamp(T_ground, T_air - 25.0, T_air + 25.0) T_snow1 = T_snow T_snow2 = T_snow elseif depth_snow > 0.05 - ttt = cp_ice * density_snow * 0.02 / length_step + ttt = cp_ice * ρ_snow * 0.02 / length_step - T_snow = (T_snow_last * ttt * ra_g * 0.04 + Gg * ra_g * 0.02 + density_air * cp_air * T_air * 0.04 + ra_g * lambda_snow * T_snow1_last) / - (density_air * cp_air * 0.04 + ra_g * lambda_snow + ttt * ra_g * 0.04) + T_snow = (T_snow_last * ttt * ra_g * 0.04 + Gg * ra_g * 0.02 + ρ_air * cp * T_air * 0.04 + ra_g * λ_snow * T_snow1_last) / + (ρ_air * cp * 0.04 + ra_g * λ_snow + ttt * ra_g * 0.04) T_snow = clamp(T_snow, T_air - 25.0, T_air + 25.0) - heat_flux_snow = lambda_snow * (T_snow - T_snow1_last) / 0.04 + G_snow = λ_snow * (T_snow - T_snow1_last) / 0.04 - heat_flux = heat_flux_snow - heat_flux = clamp(heat_flux, -100.0, 100.0) + G = G_snow + G = clamp(G, -100.0, 100.0) - heat_flux_snow1 = lambda_snow * (T_snow1_last - T_snow2_last) / (depth_snow - 0.02) - T_snow1 = T_snow1_last + ((heat_flux - heat_flux_snow1) / (cp_ice * density_snow * 0.02)) * length_step - heat_flux_snow2 = (T_snow2_last - T_any0_last) / ((0.5 * (depth_snow - 0.04) / lambda_snow) + (0.02 / lambda_soil1)) - T_snow2 = T_snow2_last + ((heat_flux_snow1 - heat_flux_snow2) / (cp_ice * density_snow * (depth_snow - 0.04))) * length_step + G_snow1 = λ_snow * (T_snow1_last - T_snow2_last) / (depth_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 - T_any0 = T_any0_last + ((heat_flux_snow2 - heat_flux_soil1) / (capacity_heat_soil0 * 0.02)) * length_step + T_any0 = T_any0_last + ((G_snow2 - G_soil1) / (capacity_heat_soil0 * 0.02)) * length_step T_soil0 = T_any0 if T_snow > zero && T_snow_last <= zero && depth_snow > zero @@ -115,6 +116,6 @@ function surface_temperature_jl(T_air::FT, rh_air::FT, depth_snow::FT, depth_wat T_ground = T_snow end - return heat_flux, T_ground, T_any0, T_soil0, T_snow, T_snow1, T_snow2 + return G, T_ground, T_any0, T_soil0, T_snow, T_snow1, T_snow2 end diff --git a/src/DataType/InterTempVars.jl b/src/DataType/InterTempVars.jl index 133b085..5a3cf25 100644 --- a/src/DataType/InterTempVars.jl +++ b/src/DataType/InterTempVars.jl @@ -1,8 +1,8 @@ @with_kw mutable struct InterTempVars Tc_u::Vector{FT} = zeros(MAX_Loop) Ts0::Vector{FT} = zeros(MAX_Loop) - Tsn0::Vector{FT} = zeros(MAX_Loop) Tsm0::Vector{FT} = zeros(MAX_Loop) + Tsn0::Vector{FT} = zeros(MAX_Loop) Tsn1::Vector{FT} = zeros(MAX_Loop) Tsn2::Vector{FT} = zeros(MAX_Loop) Qhc_o::Vector{FT} = zeros(MAX_Loop) @@ -34,6 +34,9 @@ Evap_SW::Vector{FT} = zeros(MAX_Loop) Evap_SS::Vector{FT} = zeros(MAX_Loop) lambda_snow::Vector{FT} = zeros(MAX_Loop) + + # 记录土壤温度 + Cs::Matrix{FT} = zeros(layer + 2, MAX_Loop) Tm::Matrix{FT} = zeros(layer + 2, MAX_Loop) G::Matrix{FT} = zeros(layer + 2, MAX_Loop) diff --git a/src/beps_inter_prg.jl b/src/beps_inter_prg.jl index 659e978..7cea9d0 100644 --- a/src/beps_inter_prg.jl +++ b/src/beps_inter_prg.jl @@ -17,10 +17,10 @@ The inter-module function between main program and modules """ function inter_prg_jl( jday::Int, rstep::Int, - lai::T, clumping::T, parameter::Vector{T}, meteo::ClimateData, CosZs::T, + lai::T, clumping::T, param::Vector{T}, meteo::ClimateData, CosZs::T, var_o::Vector{T}, var_n::Vector{T}, soilp::Union{Soil_c,Soil}, Ra::Radiation, - mid_res::Results, mid_ET::OutputET, var::InterTempVars; debug=false, kw...) where {T} + mid_res::Results, mid_ET::OutputET, var::InterTempVars; kw...) where {T} # var = var2 # var = InterTempVars() @@ -73,58 +73,51 @@ function inter_prg_jl( Kn = 0.3 #0.713/2.4 G_theta = 0.5 gamma = 0.066 - # double K,Vcmax0,Vcmax_sunlit,Vcmax_shaded,expr1,expr2,expr3; - # double slope_Vcmax_N, leaf_N,Jmax_sunlit,Jmax_shaded; - alpha_sat = parameter[24+1] # albedo of saturated/dry soil for module rainfall 1 - alpha_dry = parameter[25+1] # the albedo of dry soil - canopyh_o = parameter[29+1] # to be used for module aerodynamic_conductance - canopyh_u = parameter[30+1] - height_wind_sp = parameter[31+1] # the_height_to_measure_wind_speed, for module aerodynamic_conductance - m_h2o = parameter[33+1] # to be used for module photosynthesis - b_h2o = parameter[34+1] + alpha_sat = param[24+1] # albedo of saturated/dry soil for module rainfall 1 + alpha_dry = param[25+1] # the albedo of dry soil + canopyh_o = param[29+1] # to be used for module aerodynamic_conductance + canopyh_u = param[30+1] + height_wind_sp = param[31+1] # the_height_to_measure_wind_speed, for module aerodynamic_conductance + m_h2o = param[33+1] # to be used for module photosynthesis + b_h2o = param[34+1] # /***** Vcmax-Nitrogen calculations，by G.Mo，Apr. 2011 *****/ if (CosZs > 0) # day time K = G_theta * clumping / CosZs # G_theta = 0.5 assuming a spherical leaf angle distribution - Vcmax0 = parameter[36+1] + Vcmax0 = param[36+1] expr1 = 1 - exp(-K * lai) expr2 = 1 - exp(-lai * (Kn + K)) expr3 = 1 - exp(-Kn * lai) # Formulas based on Chen et al., 2012, GBC if (expr1 > 0) - Vcmax_sunlit = Vcmax0 * parameter[47+1] * parameter[46+1] * K * expr2 / (Kn + K) / expr1 + Vcmax_sunlit = Vcmax0 * param[47+1] * param[46+1] * K * expr2 / (Kn + K) / expr1 else Vcmax_sunlit = Vcmax0 end if (K > 0 && lai > expr1 / K) - Vcmax_shaded = Vcmax0 * parameter[47+1] * parameter[46+1] * (expr3 / Kn - expr2 / (Kn + K)) / (lai - expr1 / K) + Vcmax_shaded = Vcmax0 * param[47+1] * param[46+1] * (expr3 / Kn - expr2 / (Kn + K)) / (lai - expr1 / K) else Vcmax_shaded = Vcmax0 end end # /***** LAI calculation module, by B. Chen *****/ - lai_o = lai - if (lai < 0.1) - lai_o = 0.1 - end - landcover = Int(parameter[4+1]) - + lai_o = lai < 0.1 ? 0.1 : lai + + landcover = Int(param[4+1]) if (landcover == 25 || landcover == 40) lai_u = 0.01 else lai_u = 1.18 * exp(-0.99 * lai_o) end - if (lai_u > lai_o) - lai_u = 0.01 - end - - stem_o = parameter[8+1] * 0.2 # parameter[8].LAI max overstory - stem_u = parameter[9+1] * 0.2 # parameter[9].LAI max understory + (lai_u > lai_o) && (lai_u = 0.01) + + stem_o = param[8+1] * 0.2 # parameter[8].LAI max overstory + stem_u = param[9+1] * 0.2 # parameter[9].LAI max understory # lai2: separate lai into sunlit and shaded portions LAI = Leaf() @@ -138,13 +131,12 @@ function inter_prg_jl( prcp = meteo.rain / step # precipitation in meters Ta = meteo.temp - VPS_air = cal_es(Ta) # to estimate saturated water vapor pressure in kpa - e_a10 = VPS_air * rh_air / 100 # to be used for module photosynthesis - VPD_air = VPS_air - e_a10 # water vapor deficit at the reference height + es = cal_es(Ta) # to estimate saturated water vapor pressure in kpa + ea = es * rh_air / 100 # to be used for module photosynthesis + VPD_air = es - ea # water vapor deficit at the reference height - q_ca = 0.622 * e_a10 / (101.35 - 0.378 * e_a10) # in g/g, unitless + q_ca = 0.622 * ea / (101.35 - 0.378 * ea) # in g/g, unitless Cp_ca = Cpd * (1 + 0.84 * q_ca) - slope = cal_slope(Ta) if (Ks <= 0) @@ -153,54 +145,45 @@ function inter_prg_jl( alpha_v_u = 0.0 alpha_n_u = 0.0 else - alpha_v_o = parameter[22+1] - alpha_n_o = parameter[23+1] - alpha_v_u = parameter[22+1] - alpha_n_u = parameter[23+1] + alpha_v_o = param[22+1] + alpha_n_o = param[23+1] + alpha_v_u = param[22+1] + alpha_n_u = param[23+1] end - # Ground surface temperature var.Ts0[1] = clamp(var_o[3+1], Ta - 2.0, Ta + 2.0) - var.Tsn0[1] = clamp(var_o[4+1], Ta - 2.0, Ta + 2.0) var.Tsm0[1] = clamp(var_o[5+1], Ta - 2.0, Ta + 2.0) - var.Tsn1[1] = clamp(var_o[6+1], Ta - 2.0, Ta + 2.0) - var.Tsn2[1] = clamp(var_o[7+1], Ta - 2.0, Ta + 2.0) + var.Tsn0[1] = clamp(var_o[4+1], Ta - 2.0, Ta + 2.0) # snow0 + var.Tsn1[1] = clamp(var_o[6+1], Ta - 2.0, Ta + 2.0) # snow1 + var.Tsn2[1] = clamp(var_o[7+1], Ta - 2.0, Ta + 2.0) # snow2 var.Qhc_o[1] = var_o[11+1] var.Wcl_o[1] = var_o[15+1] - var.Wcs_o[1] = var_o[16+1] # /* the mass of intercepted liquid water and snow, overstory */ + var.Wcs_o[1] = var_o[16+1] # the mass of intercepted liquid water and snow, overstory # the evaporation rate of rain and snow--in kg/m^2/s, understory var.Wcl_u[1] = var_o[18+1] - var.Wcs_u[1] = var_o[19+1] # /* the mass of intercepted liquid water and snow, overstory */ - var.Wg_snow[1] = var_o[20+1] # /* thr fraction of ground surface covered in snow and snow mass */ + var.Wcs_u[1] = var_o[19+1] # the mass of intercepted liquid water and snow, overstory + var.Wg_snow[1] = var_o[20+1] # thr fraction of ground surface covered in snow and snow mass mass_water_g = Ref(0.0) Zsp = Ref(soilp.Zsp) Zp = Ref(soilp.Zp) - - if (Zp[] < 0.001) - Zp[] = 0 - end - + (Zp[] < 0.001) && (Zp[] = 0.0) + init_leaf_dbl(Tc_old, Ta - 0.5) - # Cs = zeros(layer + 2, MAX_Loop) - # Tm = zeros(layer + 2, MAX_Loop) - # G = zeros(layer + 2, MAX_Loop) - - # /***** Vcmax Jmax module by L. He *****/ + # /***** Ten time intervals in a hourly time step.6min or 360s per loop ******/ - # for(kkk = 1;kkk <= kloop;kkk++) @inbounds for kkk = 2:kloop+1 - # /***** Snow pack stage 1 by X. Luo *****/ + var.Xcs_o[kkk], var.Xcs_u[kkk], var.Xg_snow[kkk] = snowpack_stage1_jl(Ta, prcp, # var.Wcs_o[kkk-1], var.Wcs_u[kkk-1], var.Wg_snow[kkk-1], Ref(var.Wcs_o, kkk), Ref(var.Wcs_u, kkk), Ref(var.Wg_snow, kkk), lai_o, lai_u, clumping, Ref(var.Ac_snow_o, kkk), Ref(var.Ac_snow_u, kkk), Ref(var.rho_snow, kkk), Zsp, - Ref(var.alpha_v_sw, kkk), Ref(var.alpha_n_sw, kkk)) + Ref(var.alpha_v_sw, kkk), Ref(var.alpha_n_sw, kkk)) # by X. Luo # /***** Rain fall stage 1 by X. Luo *****/ var.Wcl_o[kkk], var.Wcl_u[kkk], var.Xcl_o[kkk], var.Xcl_u[kkk], var.r_rain_g[kkk] = @@ -221,7 +204,7 @@ function inter_prg_jl( # /***** Soil water factor module by L. He *****/ soil_water_factor_v2(soilp) f_soilwater = min(soilp.f_soilwater, 1.0) # used in `photosynthesis` - + GH_o = var.Qhc_o[kkk-1]# to be used as the init. for module aerodynamic_conductance init_leaf_dbl(Ci_old, 0.7 * CO2_air) @@ -268,7 +251,7 @@ function inter_prg_jl( if (CosZs > 0) photosynthesis(Tc_old, Rns, Ci_old, leleaf, - Ta, e_a10, f_soilwater, b_h2o, m_h2o, + Ta, ea, f_soilwater, b_h2o, m_h2o, Gb_o, Gb_u, Vcmax_sunlit, Vcmax_shaded, Gs_new, Ac, Ci_new; version="julia") else @@ -313,20 +296,17 @@ function inter_prg_jl( end# end of while multiply!(GPP, Ac, LAI) - # /***** Transpiration module by X. Luo *****/ - var.Trans_o[kkk], var.Trans_u[kkk] = transpiration_jl(Tc_new, Ta, rh_air, Gw, LAI) + var.Trans_o[kkk], var.Trans_u[kkk] = transpiration_jl(Tc_new, Ta, rh_air, Gw, LAI) # X. Luo # /***** Evaporation and sublimation from canopy by X. Luo *****/ var.Eil_o[kkk], var.Eil_u[kkk], var.EiS_o[kkk], var.EiS_u[kkk] = evaporation_canopy_jl( - Tc_new, Ta, rh_air, + Tc_new, Ta, rh_air, Gww, PAI, var.Xcl_o[kkk], var.Xcl_u[kkk], var.Xcs_o[kkk], var.Xcs_u[kkk]) - # /***** Rainfall stage 2 by X. Luo *****/ - rainfall_stage2_jl(var.Eil_o[kkk], var.Eil_u[kkk], Ref(var.Wcl_o, kkk), Ref(var.Wcl_u, kkk)) + rainfall_stage2_jl(var.Eil_o[kkk], var.Eil_u[kkk], Ref(var.Wcl_o, kkk), Ref(var.Wcl_u, kkk)) # X. Luo - # /***** Snow pack stage 2 by X. Luo *****/ - snowpack_stage2_jl(var.EiS_o[kkk], var.EiS_u[kkk], Ref(var.Wcs_o, kkk), Ref(var.Wcs_u, kkk)) + snowpack_stage2_jl(var.EiS_o[kkk], var.EiS_u[kkk], Ref(var.Wcs_o, kkk), Ref(var.Wcs_u, kkk)) # X. Luo # /***** Evaporation from soil module by X. Luo *****/ Gheat_g = 1 / ra_g @@ -349,15 +329,14 @@ function inter_prg_jl( var.Cs[1, kkk] = soilp.Cs[1] # added var.Cs[2, kkk] = soilp.Cs[1] - var.Tc_u[kkk] = Tcu # added - lambda[2] = soilp.lambda[1] - d_soil[2] = soilp.d_soil[1] + var.Tc_u[kkk] = Tcu # added + lambda[2] = soilp.lambda[1] + d_soil[2] = soilp.d_soil[1] var.Tm[1, kkk-1] = soilp.temp_soil_p[1] var.Tm[2, kkk-1] = soilp.temp_soil_p[2] # first place is temp_soil_p[0]? var.G[2, kkk] = soilp.G[1] - ## 二维数组`Ref`如何处理？ var.G[1, kkk], var.Ts0[kkk], var.Tm[1, kkk], var.Tsm0[kkk], var.Tsn0[kkk], var.Tsn1[kkk], var.Tsn2[kkk] = surface_temperature_jl(Ta, rh_air, Zsp[], Zp[], @@ -371,19 +350,11 @@ function inter_prg_jl( Update_temp_soil_c(soilp, var.Tm[1, kkk]) # soilp.temp_soil_c[1] = var.Tm[1, kkk] - # /***** Snow Pack Stage 3 module by X. Luo *****/ - snowpack_stage3_jl(Ta, var.Tsn0[kkk], var.Tsn0[kkk-1], var.rho_snow[kkk], Zsp, Zp, Ref(var.Wg_snow, kkk)) + snowpack_stage3_jl(Ta, var.Tsn0[kkk], var.Tsn0[kkk-1], var.rho_snow[kkk], Zsp, Zp, Ref(var.Wg_snow, kkk)) # X. Luo - # /***** Sensible heat flux module by X. Luo *****/ var.Qhc_o[kkk], var.Qhc_u[kkk], var.Qhg[kkk] = sensible_heat_jl(Tc_new, var.Ts0[kkk], Ta, rh_air, - Gh, Gheat_g, PAI) - - # println("===========================") - # println("kkk = $kkk") - # @printf("var.Ts0 = %f, var.Tm = %f, Tsno = %f, var.Tsm0 = %f, var.Tsn1 = %f, var.Tsn2 = %f, var.G = %f\n", - # var.Ts0[kkk], var.Tm[1, kkk], var.Tsn0[kkk], var.Tsm0[kkk], var.Tsn1[kkk], var.Tsn2[kkk], var.G[1, kkk]) - # @printf("Ta = %f, Gg = %f, QHs = %f, %f, %f\n", Ta, Gheat_g, var.Qhc_o[kkk], var.Qhc_u[kkk], var.Qhg[kkk]) + Gh, Gheat_g, PAI) # X. Luo # /***** Soil water module by L. He *****/ soilp.Zsp = Zsp[]