修修删删, add examples

examples/BEPS_compare_with_C.qmd

@@ -0,0 +1,61 @@
+# BEPS
+
+```{julia}
+# import RTableTools: fread
+# d = fread("examples/input/p1_meteo")
+```
+
+```{julia}
+using BEPS
+d = deserialize("data/p1_meteo")
+lai = readdlm("examples/input/p1_lai.txt")[:]
+
+par = (lon=120.5, lat=30.5, landcover=25, clumping=0.85,
+  soil_type=8, Tsoil=2.2,
+  soilwater=0.4115, snowdepth=0.0)
+```
+
+```{julia}
+@time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
+@time df_c, df_ET_c = besp_main(d, lai, par; version="c");
+
+sum(df_jl)
+sum(df_c)
+
+df_diff = abs.(df_jl .- df_c)
+df_diff_perc = abs.(df_jl .- df_c) ./ df_c .* 100
+# df_diff = abs.(df_ET_jl .- df_ET_c)
+# fwrite(df_diff, "a-3.csv")
+max(df_diff)
+max(df_diff_perc) # SH, 1.48%的误差
+```
+
+# 绘图展示结果
+```{julia}
+using Plots
+using Ipaper
+
+function plot_var(var=:SH)
+  n = size(df_jl, 1)
+  inds = 1:n
+
+  y_jl = df_jl[inds, var]
+  y_c = df_c[inds, var]
+  diff_perc = (y_jl .- y_c) ./ y_c .* 100
+
+  plot(diff_perc; title="$var", label="bias (%)") # 
+end
+
+vars = names(df_jl)[[1:4; 8; 12:16]]
+ps = map(plot_var, vars)
+plot(ps..., layout=(3, 4), size=(1400, 800))
+# write_fig("images/Figure1_bias_of_julia-version.png", 15, 8; show=false)
+```
+
+```{julia}
+# using ProfileView
+# @time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
+# @profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
+@profview_allocs df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
+@profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
+```

examples/example_01.qmd

@@ -24,127 +24,32 @@ sum(df_c)
 
 df_diff = abs.(df_jl .- df_c)
 df_diff_perc = abs.(df_jl .- df_c) ./ df_c .* 100
-# df_diff = abs.(df_ET_jl .- df_ET_c)
-# fwrite(df_diff, "a-3.csv")
+
 max(df_diff)
 max(df_diff_perc) # SH, 1.48%的误差
 ```
 
 # 绘图展示结果
+
 ```{julia}
 using Plots
 using Ipaper
+gr(framestyle=:box)
 
 function plot_var(var=:SH)
   n = size(df_jl, 1)
+  n = 24*7
   inds = 1:n
 
   y_jl = df_jl[inds, var]
-  y_c = df_c[inds, var]
-  diff_perc = (y_jl .- y_c) ./ y_c .* 100
-
-  plot(diff_perc; title="$var", label="bias (%)") # 
+  # y_c = df_c[inds, var]
+  # diff_perc = (y_jl .- y_c) ./ y_c .* 100
+  plot(y_jl; title="$var", label=nothing) # 
 end
 
 vars = names(df_jl)[[1:4; 8; 12:16]]
 ps = map(plot_var, vars)
-# p1 = plot(df_jl.SH[inds], label="julia", title="SH")
-# plot!(df_c.SH[inds] + 0, label="c")
-# p2 = plot(df_diff_perc.SH[inds], label="Julia error (%)")
-# title = "bias (%): (julia - c)/c"
-plot(ps..., layout=(3, 4), size=(1400, 800))
-write_fig("images/Figure1_bias_of_julia-version.png", 15, 8; show=false)
-```
-
-```{julia}
-# using ProfileView
-# @time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
-# @profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
-@profview_allocs df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
-@profview df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
-```
-
-
-```{julia}
-## C版本的结果
-# :gpp_o_sunlit => 41704.09054753295
-# :gpp_u_sunlit => 285.1932406544977
-# :gpp_o_shaded => 12720.54113242912
-# :gpp_u_shaded => 135.13475998737647
-# :plant_resp => 0.0
-# :npp_o => 0.0
-# :npp_u => 0.0
-# :GPP => 2369.3022582020903
-# :NPP => 0.0
-# :NEP => 0.0
-# :soil_resp => 0.0
-# :Net_Rad => 1.299620168997007e6
-# :SH => 232161.16420363513
-# :LH => 811896.0668881282
-# :Trans => 444.02437414775443
-# :Evap => 748.8880942489606
 
-## sum
-:gpp_o_sunlit => 0.0002343150940907193
-:gpp_u_sunlit => 0.05115824940432174
-:gpp_o_shaded => 1.4900267617763002e-5
-:gpp_u_shaded => 1.0758724368703184e-7
-:plant_resp   => 0.0
-:npp_o        => 0.0
-:npp_u        => 0.0
-:GPP          => 0.0022203403965193147
-:NPP          => 0.0
-:NEP          => 0.0
-:soil_resp    => 0.0
-:Net_Rad      => 0.003708570688994328
-:SH           => 0.12619010825226695
-:LH           => 2.5288028670446403
-:Trans        => 0.0035486520015893166
-:Evap         => 0.0009524622286666777
-
-## max
- :gpp_o_sunlit => 0.00021630367300851105
- :gpp_u_sunlit => 0.0005586954093889651
- :gpp_o_shaded => 9.06128597910616e-6
- :gpp_u_shaded => 1.7276561885964936e-8
-   :plant_resp => 0.0
-        :npp_o => 0.0
-        :npp_u => 0.0
-          :GPP => 2.4137032993143404e-5
-          :NPP => 0.0
-          :NEP => 0.0
-    :soil_resp => 0.0
-      :Net_Rad => 0.002326329983901587
-           :SH => 0.029240292020906722
-           :LH => 0.07646965953554741
-        :Trans => 1.6816338803168907e-5
-         :Evap => 0.00012368929687477104
-```
-
-
-# Test about
-
-```{julia}
-canopyh_o = 2.0
-canopyh_u = 0.2
-height_wind_sp = 2.0
-clumping = 0.8
-Ta = 20.0
-wind_sp = 2.
-GH_o = 100.0
-pai_o = 4.0
-pai_u = 2.0
-
-ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
-  aerodynamic_conductance_c(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
-    pai_o, pai_u)
-r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)
-
-ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
-  aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
-    pai_o, pai_u)
-r2 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)
-
-@show r1
-@show r2;
+plot(ps..., layout=(3, 4), size=(1400, 800))
+# write_fig("images/Figure1_bias_of_julia-version.png", 15, 8; show=false)
 ```

src/BEPS/aerodynamic_conductance.jl

@@ -5,8 +5,20 @@
 
 # Examples
 ```julia
-ra_o, G_o_a, G_o_b, G_u_a, G_u_b, ra_g = aerodynamic_conductance_jl(
-    canopy_height_o, canopy_height_u, z_wind, clumping, temp_air, wind_sp, SH_o_p, lai_o, lai_u)
+canopyh_o = 2.0
+canopyh_u = 0.2
+height_wind_sp = 2.0
+clumping = 0.8
+Ta = 20.0
+wind_sp = 2.
+GH_o = 100.0
+lai_o = 4.0
+lai_u = 2.0
+
+ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
+  aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
+    lai_o, lai_u)
+r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)
 ```
 """
 function aerodynamic_conductance_jl(canopy_height_o::FT, canopy_height_u::FT,
@@ -91,7 +103,7 @@ end
 
 
 function windProfile_factor(canopy_height_u, canopy_height_o, gamma, k=1.0)
-  exp(-gamma * (1 - canopy_height_u * k / canopy_height_o))
+  exp(gamma * (1 - canopy_height_u * k / canopy_height_o))
 end
 
 function cal_Nu(u::FT, nu_lower::FT)::FT

src/beps_inter_prg.jl

@@ -36,10 +36,12 @@ function inter_prg_jl(
   Tc_old = Leaf()
   Tc_new = Leaf()
   Gs_old = Leaf()
-  Gc = Leaf() # total conductance for CO2 from the intercellular space of the leaves to the reference height above the canopy
-  Gh = Leaf() # total conductance for heat transfer from the leaf surface to the reference height above the canopy
-  Gw = Leaf() # total conductance for water from the intercellular space of the leaves to the reference height above the canopy
-  Gww = Leaf() # total conductance for water from the surface of the leaves to the reference height above the canopy
+
+  # to the reference height above the canopy
+  Gc = Leaf()  # total conductance for CO2 from the intercellular space of the leaves
+  Gh = Leaf()  # total conductance for heat transfer from the leaf surface 
+  Gw = Leaf()  # total conductance for water from the intercellular space of the leaves
+  Gww = Leaf() # total conductance for water from the surface of the leaves
 
   Gs_new = Leaf()
   Ac = Leaf()
@@ -71,7 +73,6 @@ 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;
 
@@ -435,6 +436,5 @@ function inter_prg_jl(
 
   # total GPP . gC/m2/step
   mid_res.GPP = (GPP.o_sunlit + GPP.o_shaded + GPP.u_sunlit + GPP.u_shaded) * 12 * step * 0.000001
-  # return
   nothing
 end 

test/runtests.jl

@@ -4,7 +4,7 @@ using BEPS: path_proj
 
 
 
-include("test-readparam.jl")
+include("test-clang.jl")
 include("test-soil.jl")
 
 

test/test-clang.jl

@@ -0,0 +1,46 @@
+using Test
+using BEPS
+
+@testset "readparam" begin
+  lcs = [1, 2, 6, 9, 13, 40, -1]
+  for lc = lcs
+    par1 = readparam(lc)
+    par2 = clang.readparam(lc)
+    @test maximum(abs.(par1 .- par2)) == 0
+  end
+end
+
+@testset "readcoef" begin
+  lcs = [1, 6, 13, -1]
+  stxts = 1:11
+
+  for lc = lcs, stxt = stxts
+    coef1 = readcoef(lc, stxt)
+    coef2 = clang.readcoef(lc, stxt)
+    @test maximum(abs.(coef1 .- coef2)) == 0
+  end
+end
+
+@testset "aerodynamic_conductance" begin
+  canopyh_o = 2.0
+  canopyh_u = 0.2
+  height_wind_sp = 2.0
+  clumping = 0.8
+  Ta = 20.0
+  wind_sp = 2.0
+  GH_o = 100.0
+  pai_o = 4.0
+  pai_u = 2.0
+
+  ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
+    aerodynamic_conductance_jl(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
+      pai_o, pai_u)
+  r1 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)
+
+  ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u =
+    clang.aerodynamic_conductance_c(canopyh_o, canopyh_u, height_wind_sp, clumping, Ta, wind_sp, GH_o,
+      pai_o, pai_u)
+  r2 = (; ra_o, ra_u, ra_g, Ga_o, Gb_o, Ga_u, Gb_u)
+
+  @test maximum(abs.(values(r1) .- values(r2))) <= 1e-8
+end

test/test-readparam.jl

@@ -2,24 +2,3 @@ using Test
 using BEPS
 # import BEPS: readparam, readcoef
 
-@testset "readparam" begin
-  lcs = [1, 2, 6, 9, 13, 40, -1]
-
-  for lc = lcs
-    par1 = readparam(lc)
-    par2 = clang.readparam(lc)
-    @test maximum(abs.(par1 .- par2)) == 0
-  end
-end
-
-
-@testset "readparam" begin
-  lcs = [1, 6, 13, -1]
-  stxts = 1:11
-
-  for lc = lcs, stxt = stxts
-    coef1 = readcoef(lc, stxt)
-    coef2 = clang.readcoef(lc, stxt)
-    @test maximum(abs.(coef1 .- coef2)) == 0
-  end
-end