separate cal_Rln_Longwave

jl-pkgs · Jan 16, 2024 · db6c78c · db6c78c
1 parent e78b555
commit db6c78c
Show file tree

Hide file tree

Showing 7 changed files with 148 additions and 163 deletions.
diff --git a/scripts/s1_GEE.qmd → data/s1_GEE.qmd b/scripts/s1_GEE.qmd → data/s1_GEE.qmd
diff --git 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
@@ -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
@@ -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
@@ -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,29 +12,29 @@ 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
   T_soil0::FT = 0.0
   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
@@ -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)