diff --git a/src/BEPS.jl b/src/BEPS.jl
index fd77dad..8f1408c 100644
--- a/src/BEPS.jl
+++ b/src/BEPS.jl
@@ -5,10 +5,12 @@ using UnPack
 import Parameters: @with_kw, @with_kw_noshow
 using Printf
 using Reexport
+using DataFrames: DataFrame
+
 @reexport using Serialization: deserialize, serialize
 @reexport using DelimitedFiles: readdlm
 export besp_main
-export init_soil!
+export Init_Soil_var_o
 export Soil_c
 
 path_proj(f...) = normpath(joinpath(@__DIR__, "..", f...))
@@ -26,8 +28,6 @@ include("clang/BEPS_c.jl")
 @reexport import BEPS.clang;
 import BEPS.clang: inter_prg_c, photosynthesis_c, Soil_c
 
-include("BEPS_helper.jl")
-
 include("beps_inter_prg.jl")
 include("beps_main.jl")
 
diff --git a/src/BEPS/Base/Base.jl b/src/BEPS/Base/Base.jl
index 0c37dcc..3e8c77c 100644
--- a/src/BEPS/Base/Base.jl
+++ b/src/BEPS/Base/Base.jl
@@ -1,11 +1,14 @@
 using DocStringExtensions: TYPEDFIELDS
+using DataFrames: DataFrame
 
 include("Ipaper.jl")
 include("helpers.jl")
 include("s_coszs.jl")
-include("lai2.jl")
-include("aerodynamic_conductance.jl")
 
 include("AirLayer.jl")
+include("fill_meteo.jl")
+
+include("lai2.jl")
+include("aerodynamic_conductance.jl")
 
-export s_coszs, lai2, aerodynamic_conductance_jl
+export s_coszs, lai2, aerodynamic_conductance_jl, fill_meteo!
diff --git a/src/BEPS/Base/fill_meteo.jl b/src/BEPS/Base/fill_meteo.jl
new file mode 100644
index 0000000..59f0757
--- /dev/null
+++ b/src/BEPS/Base/fill_meteo.jl
@@ -0,0 +1,31 @@
+using DataFrames: DataFrame
+
+
+function q2RH(hum::FT, tem::FT)::FT
+  # Vapour pressure in mbar
+  ea = 0.46 * hum * (tem + 273.16) / 100
+  es = 6.1078 * exp((17.269 * tem) / (237.3 + tem))
+  clamp(ea / es * 100, 0.0, 100.0)
+end
+
+"""
+Srad: W m-2
+temp: degC
+rain: m
+wind: m/s
+hum: specific humidity, q
+"""
+function fill_meteo!(meteo::ClimateData,
+  rad::FT, tem::FT, pre::FT, wind::FT, hum::FT)
+
+  meteo.Srad = rad
+  meteo.temp = tem
+  meteo.rain = pre / 1000 # m to mm
+  meteo.wind = wind
+  meteo.LR = -200.0 #  -200.0 means no measured long-wave radiation, the value will be 
+  meteo.rh = q2RH(hum, tem)
+end
+
+function fill_meteo!(meteo::ClimateData, d::DataFrame, k::Int=1)
+  fill_meteo!(meteo, d.rad[k], d.tem[k], d.pre[k], d.wind[k], d.hum[k])
+end
diff --git a/src/BEPS/Soil/Init_Soil_Parameters.jl b/src/BEPS/Soil/Init_Soil_Parameters.jl
index f75eb9f..6eb8a75 100644
--- a/src/BEPS/Soil/Init_Soil_Parameters.jl
+++ b/src/BEPS/Soil/Init_Soil_Parameters.jl
@@ -1,3 +1,6 @@
+
+
+
 function Init_Soil_Parameters(landcover::Integer, stxt::Integer, r_root_decay::Float64, p::Soil)
   p.n_layer = 5
 
diff --git a/src/BEPS/Soil/Init_Soil_var_o.jl b/src/BEPS/Soil/Init_Soil_var_o.jl
new file mode 100644
index 0000000..2aefe60
--- /dev/null
+++ b/src/BEPS/Soil/Init_Soil_var_o.jl
@@ -0,0 +1,30 @@
+function Init_Soil_var_o(p_soil::AbstractSoil, var_o::Vector,
+  meteo::ClimateData, parameter::Vector, par::NamedTuple)
+  # landcover::Int, soil_type::Int, Tsoil, soilwater, snowdepth
+  # /***** initialize soil conditions, read soil parameters and set depth *****/
+  Init_Soil_Parameters(par.landcover, par.soil_type, parameter[28], p_soil)
+  p_soil.r_drainage = parameter[27]
+  Init_Soil_Status(p_soil, par.Tsoil, meteo.temp, par.soilwater, par.snowdepth) # LHE
+
+  # Initialize intermediate variables array
+  var_o .= 0
+  for i = 4:9
+    var_o[i] = meteo.temp
+  end
+  for i = 10:15
+    var_o[i] = p_soil.Tsoil_p[i-9]
+  end
+  for i = 22:27
+    var_o[i] = p_soil.θ_prev[i-21]
+  end
+  for i = 28:33
+    var_o[i] = p_soil.ice_ratio[i-27]
+  end
+  # for (i=0;i<=40;i++)   var_o[i] = 0;
+  # for (i=3;i<=8;i++)   var_o[i] = tem;
+  # for(i=9;i<=14;i++) var_o[i] = p_soil->Tsoil_p[i-9];
+  # for(i=21;i<=26;i++) var_o[i] = p_soil->θ_prev[i-21];
+  # for(i=27;i<=32;i++) var_o[i] = p_soil->ice_ratio[i-27];
+  nothing
+end
+
diff --git a/src/BEPS/Soil/Soil.jl b/src/BEPS/Soil/Soil.jl
index 51e4faf..8e3c2b8 100644
--- a/src/BEPS/Soil/Soil.jl
+++ b/src/BEPS/Soil/Soil.jl
@@ -1,4 +1,5 @@
 include("Init_Soil_Parameters.jl")
+include("Init_Soil_var_o.jl")
 include("Update_thermal.jl")
 include("UpdateSoilMoisture.jl")
 
diff --git a/src/BEPS/snowpack.jl b/src/BEPS/snowpack.jl
index 8fecc9a..313a908 100644
--- a/src/BEPS/snowpack.jl
+++ b/src/BEPS/snowpack.jl
@@ -3,7 +3,7 @@ function snowpack_stage1_jl(temp_air::Float64, prcp::Float64,
   mass_snow_o::Ref{Float64}, mass_snow_u::Ref{Float64}, mass_snow_g::Ref{Float64},
   lai_o::Float64, lai_u::Float64, clumping::Float64, area_snow_o::Ref{Float64}, area_snow_u::Ref{Float64},
   # perc_snow_o::Ref{Float64}, perc_snow_u::Ref{Float64}, perc_snow_g::Ref{Float64},
-  density_snow::Ref{Float64}, depth_snow::Ref{Float64}, albedo_v_snow::Ref{Float64}, albedo_n_snow::Ref{Float64})
+  ρ_snow::Ref{Float64}, depth_snow::Ref{Float64}, albedo_v_snow::Ref{Float64}, albedo_n_snow::Ref{Float64})
 
   massMax_snow_o = 0.1 * lai_o
   massMax_snow_u = 0.1 * lai_u
@@ -12,8 +12,8 @@ function snowpack_stage1_jl(temp_air::Float64, prcp::Float64,
 
   length_step = kstep
 
-  density_new_snow = 67.9 + 51.3 * exp(temp_air / 2.6)
-  density_water = 1025
+  ρ_new_snow = 67.9 + 51.3 * exp(temp_air / 2.6)
+  ρ_water = 1025
   albedo_v_Newsnow = 0.94
   albedo_n_Newsnow = 0.8
 
@@ -21,25 +21,25 @@ function snowpack_stage1_jl(temp_air::Float64, prcp::Float64,
   mass_snow_u_last = mass_snow_u[]
   mass_snow_g_last = mass_snow_g[]
 
-  snowrate = temp_air > 0 ? 0 : prcp * density_water / density_new_snow
+  snowrate = temp_air > 0 ? 0 : prcp * ρ_water / ρ_new_snow
   snowrate_o = 0.0
 
   if temp_air < 0
     snowrate_o = snowrate
-    mass_snow_o[] = mass_snow_o_last + snowrate_o * length_step * density_new_snow * (1 - exp(-lai_o * clumping))
+    mass_snow_o[] = mass_snow_o_last + snowrate_o * length_step * ρ_new_snow * (1 - exp(-lai_o * clumping))
     perc_snow_o = mass_snow_o[] / massMax_snow_o
     perc_snow_o = clamp(perc_snow_o, 0, 1)
     area_snow_o[] = perc_snow_o * areaMax_snow_o
     massStep_snow_o = mass_snow_o[] - mass_snow_o_last
 
-    snowrate_u = max(0, snowrate_o - (massStep_snow_o) / density_new_snow / length_step)
-    mass_snow_u[] = mass_snow_u_last + snowrate_u * length_step * density_new_snow * (1 - exp(-lai_u * clumping))
+    snowrate_u = max(0, snowrate_o - (massStep_snow_o) / ρ_new_snow / length_step)
+    mass_snow_u[] = mass_snow_u_last + snowrate_u * length_step * ρ_new_snow * (1 - exp(-lai_u * clumping))
     perc_snow_u = mass_snow_u[] / massMax_snow_u
     perc_snow_u = clamp(perc_snow_u, 0, 1)
     area_snow_u[] = perc_snow_u * areaMax_snow_u
     massStep_snow_u = mass_snow_u[] - mass_snow_u_last
 
-    snowrate_g = max(0, snowrate_u - (massStep_snow_u) / density_new_snow / length_step)
+    snowrate_g = max(0, snowrate_u - (massStep_snow_u) / ρ_new_snow / length_step)
     change_depth_snow = snowrate_g * length_step
   else
     mass_snow_o[] = mass_snow_o_last
@@ -55,17 +55,17 @@ function snowpack_stage1_jl(temp_air::Float64, prcp::Float64,
   end
 
   change_depth_snow = max(0, change_depth_snow)
-  mass_snow_g[] = mass_snow_g_last + change_depth_snow * density_new_snow
+  mass_snow_g[] = mass_snow_g_last + change_depth_snow * ρ_new_snow
   mass_snow_g[] = max(0, mass_snow_g[])
 
   if change_depth_snow > 0
-    density_snow[] = (density_snow[] * depth_snow[] + density_new_snow * change_depth_snow) / (depth_snow[] + change_depth_snow)
+    ρ_snow[] = (ρ_snow[] * depth_snow[] + ρ_new_snow * change_depth_snow) / (depth_snow[] + change_depth_snow)
   else
-    density_snow[] = (density_snow[] - 250) * exp(-0.001 * length_step / 3600.0) + 250.0
+    ρ_snow[] = (ρ_snow[] - 250) * exp(-0.001 * length_step / 3600.0) + 250.0
   end
 
-  depth_snow[] = mass_snow_g[] > 0 ? mass_snow_g[] / density_snow[] : 0.0
-  perc_snow_g = min(mass_snow_g[] / (0.05 * density_snow[]), 1.0)
+  depth_snow[] = mass_snow_g[] > 0 ? mass_snow_g[] / ρ_snow[] : 0.0
+  perc_snow_g = min(mass_snow_g[] / (0.05 * ρ_snow[]), 1.0)
 
   if snowrate_o > 0
     albedo_v_snow[] = (albedo_v_snow[] - 0.70) * exp(-0.005 * length_step / 3600) + 0.7
@@ -89,7 +89,7 @@ function snowpack_stage2_jl(evapo_snow_o::Float64, evapo_snow_u::Float64,
 end
 
 
-function snowpack_stage3_jl(temp_air::Float64, temp_snow::Float64, temp_snow_last::Float64, density_snow::Float64,
+function snowpack_stage3_jl(temp_air::Float64, temp_snow::Float64, temp_snow_last::Float64, ρ_snow::Float64,
   depth_snow::Ref{Float64}, depth_water::Ref{Float64}, mass_snow_g::Ref{Float64})
 
   length_step = kstep  # length_step in model
@@ -100,8 +100,8 @@ function snowpack_stage3_jl(temp_air::Float64, temp_snow::Float64, temp_snow_las
 
   # parameters
   cp_ice = 2228.261
-  latent_fusion = 3.34 * 1000000
-  density_water = 1025
+  latent_fusion = 3.34 * 1000000 # J Kg-1
+  ρ_water = 1025 # Kg m-3
 
   # simulate snow melt and freeze process
   mass_snow_melt = 0.0
@@ -118,30 +118,30 @@ function snowpack_stage3_jl(temp_air::Float64, temp_snow::Float64, temp_snow_las
     # case 2 depth of snow > 0.02 < 0.05 m
   elseif depth_snow_sup > 0.02 && depth_snow_sup <= 0.05
     if temp_snow > 0 && temp_snow_last < 0 && mass_snow_sup > 0  # melting
-      mass_snow_melt = temp_snow * cp_ice * density_snow * depth_snow_sup / latent_fusion
+      mass_snow_melt = temp_snow * cp_ice * ρ_snow * depth_snow_sup / latent_fusion
       mass_snow_melt = min(mass_snow_sup, mass_snow_melt)  # the amount of melted snow could not be larger than supply
     else
       mass_snow_melt = 0
     end
 
     if temp_snow <= 0 && temp_snow_last > 0 && depth_water[] > 0  # freezing
-      mass_water_frozen = (0 - temp_snow) * cp_ice * density_snow * depth_snow_sup / latent_fusion
-      mass_water_frozen = max(mass_water_frozen, depth_water[] * density_water)
+      mass_water_frozen = (0 - temp_snow) * cp_ice * ρ_snow * depth_snow_sup / latent_fusion
+      mass_water_frozen = max(mass_water_frozen, depth_water[] * ρ_water)
     else
       mass_water_frozen = 0
     end
     # case 3 depth of snow > 0.05 m
   elseif depth_snow_sup > 0.05
     if temp_snow > 0 && temp_snow_last <= 0 && mass_snow_sup > 0  # melting
-      mass_snow_melt = temp_snow * cp_ice * density_snow * depth_snow_sup * 0.02 / latent_fusion
+      mass_snow_melt = temp_snow * cp_ice * ρ_snow * depth_snow_sup * 0.02 / latent_fusion
       mass_snow_melt = min(mass_snow_sup, mass_snow_melt)  # the amount of melted snow could not be larger than supply
     else
       mass_snow_melt = 0
     end
 
     if temp_snow <= 0 && temp_snow_last > 0 && depth_water[] > 0  # freezing
-      mass_water_frozen = (0 - temp_snow) * cp_ice * density_snow * depth_snow_sup * 0.02 / latent_fusion
-      mass_water_frozen = max(mass_water_frozen, depth_water[] * density_water)
+      mass_water_frozen = (0 - temp_snow) * cp_ice * ρ_snow * depth_snow_sup * 0.02 / latent_fusion
+      mass_water_frozen = max(mass_water_frozen, depth_water[] * ρ_water)
     else
       mass_water_frozen = 0
     end
@@ -152,14 +152,14 @@ function snowpack_stage3_jl(temp_air::Float64, temp_snow::Float64, temp_snow_las
   mass_snow_g[] = max(0, mass_snow_g[])
 
   # change of depth in snow
-  melt_depth_snow = mass_snow_melt / density_snow
-  frozen_depth_snow = mass_water_frozen / density_snow
+  melt_depth_snow = mass_snow_melt / ρ_snow
+  frozen_depth_snow = mass_water_frozen / ρ_snow
   depth_snow[] = depth_snow_sup - melt_depth_snow + frozen_depth_snow
   depth_snow[] = max(0, depth_snow[])
 
   # change of depth in water
-  melt_depth_water = mass_snow_melt / density_water
-  frozen_depth_water = mass_water_frozen / density_water
+  melt_depth_water = mass_snow_melt / ρ_water
+  frozen_depth_water = mass_water_frozen / ρ_water
   depth_water[] = depth_water[] + melt_depth_water - frozen_depth_water
   depth_water[] = max(0, depth_water[])
 end
diff --git a/src/BEPS_helper.jl b/src/BEPS_helper.jl
deleted file mode 100644
index 7dc93da..0000000
--- a/src/BEPS_helper.jl
+++ /dev/null
@@ -1,82 +0,0 @@
-using DataFrames: DataFrame
-
-
-function q2RH(hum::FT, tem::FT)::FT
-  # Vapour pressure in mbar
-  ea = 0.46 * hum * (tem + 273.16) / 100
-  es = 6.1078 * exp((17.269 * tem) / (237.3 + tem))
-  clamp(ea / es * 100, 0.0, 100.0)
-end
-
-"""
-Srad: W m-2
-temp: degC
-rain: m
-wind: m/s
-hum: specific humidity, q
-"""
-function fill_meteo!(meteo::ClimateData,
-  rad::FT, tem::FT, pre::FT, wind::FT, hum::FT)
-
-  meteo.Srad = rad
-  meteo.temp = tem
-  meteo.rain = pre / 1000 # m to mm
-  meteo.wind = wind
-  meteo.LR = -200.0 #  -200.0 means no measured long-wave radiation, the value will be 
-  meteo.rh = q2RH(hum, tem)
-end
-
-function fill_meteo!(meteo::ClimateData, d::DataFrame, k::Int=1)
-  fill_meteo!(meteo, d.rad[k], d.tem[k], d.pre[k], d.wind[k], d.hum[k])
-end
-
-
-function init_soil!(p_soil::Union{Soil_c,Soil}, var_o::Vector,
-  meteo::ClimateData, parameter::Vector, par::NamedTuple)
-  # landcover::Int, soil_type::Int, Tsoil, soilwater, snowdepth
-  # /***** initialize soil conditions, read soil parameters and set depth *****/
-  Init_Soil_Parameters(par.landcover, par.soil_type, parameter[28], p_soil)
-  p_soil.r_drainage = parameter[27]
-  Init_Soil_Status(p_soil, par.Tsoil, meteo.temp, par.soilwater, par.snowdepth) # LHE
-  
-  # Initialize intermediate variables array
-  var_o .= 0
-  for i = 4:9
-    var_o[i] = meteo.temp
-  end
-  for i = 10:15
-    var_o[i] = p_soil.Tsoil_p[i-9]
-  end
-  for i = 22:27
-    var_o[i] = p_soil.θ_prev[i-21]
-  end
-  for i = 28:33
-    var_o[i] = p_soil.ice_ratio[i-27]
-  end
-  # for (i=0;i<=40;i++)   var_o[i] = 0;
-  # for (i=3;i<=8;i++)   var_o[i] = tem;
-  # for(i=9;i<=14;i++) var_o[i] = p_soil->Tsoil_p[i-9];
-  # for(i=21;i<=26;i++) var_o[i] = p_soil->θ_prev[i-21];
-  # for(i=27;i<=32;i++) var_o[i] = p_soil->ice_ratio[i-27];
-  nothing
-end
-
-
-## fill values
-const TypeDF = Union{Results,ClimateData,OutputET}
-
-## put struct into a data.frame
-function Base.getindex(x::T, i::Int)::FT where {T<:TypeDF}
-  # key = fieldnames(T)[i]
-  getfield(x, i)
-end
-
-Base.length(x::T) where {T<:TypeDF} = fieldcount(T)
-
-function fill_res!(df::DataFrame, Res::T, k::Int) where {T<:TypeDF}
-  n = length(Res)
-  for i in 1:n
-    df[k, i] = Res[i]
-  end
-  nothing
-end
diff --git a/src/DataType/DataType.jl b/src/DataType/DataType.jl
index 57d9caf..9edcc92 100644
--- a/src/DataType/DataType.jl
+++ b/src/DataType/DataType.jl
@@ -36,9 +36,29 @@ include("OUTPUT.jl")
 # export TSoil
 
 
+## fill values
+const TypeDF = Union{Results,ClimateData,OutputET}
+
+## put struct into a data.frame
+function Base.getindex(x::T, i::Int)::FT where {T<:TypeDF}
+  # key = fieldnames(T)[i]
+  getfield(x, i)
+end
+
+Base.length(x::T) where {T<:TypeDF} = fieldcount(T)
+
+function fill_res!(df::DataFrame, Res::T, k::Int) where {T<:TypeDF}
+  n = length(Res)
+  for i in 1:n
+    df[k, i] = Res[i]
+  end
+  nothing
+end
+
+
 export 
   Leaf, 
-  Soil, 
+  Soil, AbstractSoil, 
   ClimateData, Results, Cpools,
   InterTempVars,
   OutputET, Radiation
diff --git a/src/DataType/Soil.jl b/src/DataType/Soil.jl
index 2b4b9c6..f29e976 100644
--- a/src/DataType/Soil.jl
+++ b/src/DataType/Soil.jl
@@ -1,4 +1,6 @@
-@with_kw mutable struct Soil
+abstract type AbstractSoil end
+
+@with_kw mutable struct Soil <: AbstractSoil
   flag        ::Cint    = Cint(0)
   n_layer     ::Cint    = Cint(5)
   step_period ::Cint    = Cint(1)
diff --git a/src/beps_main.jl b/src/beps_main.jl
index 00454b4..85c76fd 100644
--- a/src/beps_main.jl
+++ b/src/beps_main.jl
@@ -45,7 +45,7 @@ function besp_main(d::DataFrame, lai::Vector, par::NamedTuple;
       fill_meteo!(meteo, d, k)
 
       if (flag == 0)
-        init_soil!(p_soil, var_o, meteo, param, par) # update p_soil and var_o
+        Init_Soil_var_o(p_soil, var_o, meteo, param, par) # update p_soil and var_o
       else
         var_o .= v2last
       end
diff --git a/src/clang/struct_SOIL.jl b/src/clang/struct_SOIL.jl
index ec316f9..3fc3b19 100644
--- a/src/clang/struct_SOIL.jl
+++ b/src/clang/struct_SOIL.jl
@@ -1,7 +1,7 @@
 nzero(n) = tuple(zeros(n)...) # n double zero
 const NT10 = NTuple{10,Cdouble}
 
-@with_kw mutable struct Soil_c
+@with_kw mutable struct Soil_c <: AbstractSoil
   flag        ::Cint = Cint(0)
   n_layer     ::Cint = Cint(5)
   step_period ::Cint = Cint(1)
diff --git a/test/test-soil.jl b/test/test-soil.jl
index db658d1..12c2b91 100644
--- a/test/test-soil.jl
+++ b/test/test-soil.jl
@@ -39,7 +39,7 @@ end
 end
 
 
-@testset "init_soil!" begin
+@testset "Init_Soil_var_o" begin
   p_jl = Soil()
   p_c = Soil_c()
   var_jl = zeros(41)
@@ -58,8 +58,8 @@ end
 
   param = readparam(par.landcover)      # n = 48
 
-  init_soil!(p_jl, var_jl, meteo, param, par)
-  init_soil!(p_c, var_c, meteo, param, par)
+  Init_Soil_var_o(p_jl, var_jl, meteo, param, par)
+  Init_Soil_var_o(p_c, var_c, meteo, param, par)
 
   
   funs = [Update_Cs,