diff --git a/src/BEPS/BEPS.jl b/src/BEPS/BEPS.jl
index c3bb972..8ace135 100644
--- a/src/BEPS/BEPS.jl
+++ b/src/BEPS/BEPS.jl
@@ -14,10 +14,14 @@ include("transpiration.jl")
 include("evaporation_soil.jl")
 include("evaporation_canopy.jl")
 
+include("rainfall_stage.jl")
+
 include("netRadiation.jl")
 include("photosynthesis.jl")
 
-export s_coszs, lai2, readparam, readcoef
+export s_coszs, lai2, readparam, readcoef, 
+  rainfall_stage1_jl, rainfall_stage2_jl
+
 export aerodynamic_conductance_jl, 
   sensible_heat_jl,
   latent_heat!,
diff --git a/src/BEPS/latent_heat.jl b/src/BEPS/latent_heat.jl
index 29152a3..61fe376 100644
--- a/src/BEPS/latent_heat.jl
+++ b/src/BEPS/latent_heat.jl
@@ -1,6 +1,6 @@
-function latent_heat!(leleaf::Leaf, Gw::Leaf, VPD_air, slope, Tc_old::Leaf, temp_air, rho_a, Cp_ca, psychrometer)
-  leleaf.o_sunlit = Gw.o_sunlit * (VPD_air + slope * (Tc_old.o_sunlit - temp_air)) * rho_a * Cp_ca / psychrometer
-  leleaf.o_shaded = Gw.o_shaded * (VPD_air + slope * (Tc_old.o_shaded - temp_air)) * rho_a * Cp_ca / psychrometer
-  leleaf.u_sunlit = Gw.u_sunlit * (VPD_air + slope * (Tc_old.u_sunlit - temp_air)) * rho_a * Cp_ca / psychrometer
-  leleaf.u_shaded = Gw.u_shaded * (VPD_air + slope * (Tc_old.u_shaded - temp_air)) * rho_a * Cp_ca / psychrometer
+function latent_heat!(leleaf::Leaf, Gw::Leaf, VPD, slope, Tc_old::Leaf, Tair, rho_a, Cp_ca, gamma)
+  leleaf.o_sunlit = Gw.o_sunlit * (VPD + slope * (Tc_old.o_sunlit - Tair)) * rho_a * Cp_ca / gamma
+  leleaf.o_shaded = Gw.o_shaded * (VPD + slope * (Tc_old.o_shaded - Tair)) * rho_a * Cp_ca / gamma
+  leleaf.u_sunlit = Gw.u_sunlit * (VPD + slope * (Tc_old.u_sunlit - Tair)) * rho_a * Cp_ca / gamma
+  leleaf.u_shaded = Gw.u_shaded * (VPD + slope * (Tc_old.u_shaded - Tair)) * rho_a * Cp_ca / gamma
 end
diff --git a/src/BEPS/rainfall_stage.jl b/src/BEPS/rainfall_stage.jl
new file mode 100644
index 0000000..30b8033
--- /dev/null
+++ b/src/BEPS/rainfall_stage.jl
@@ -0,0 +1,51 @@
+function rainfall_stage1_jl(Tair::Float64, prcp::Float64, mass_water_o_last::Float64, mass_water_u_last::Float64,
+  lai_o::Float64, lai_u::Float64, clumping::Float64)
+
+  length_step = kstep
+
+  # Ta > 0, otherwise it is snow fall
+  Tair <= 0.0 && (prcp = 0.0)
+
+  density_water = 1025.0
+  prcp_g = 0.0
+
+  # overstorey
+  prcp_o = prcp
+  mass_water_o = mass_water_o_last + prcp_o * length_step * density_water * (1 - exp(-lai_o * clumping))
+  massMax_water_o = 0.1 * lai_o
+
+  mass_water_o = clamp(mass_water_o, 0, massMax_water_o)
+
+  massStep_water_o = mass_water_o - mass_water_o_last
+  massStep_water_o = max(0.0, massStep_water_o)
+
+  percent_water_o = mass_water_o / massMax_water_o
+  percent_water_o = min(1.0, percent_water_o)
+
+  # understorey
+  prcp_u = prcp_o - massStep_water_o / density_water / length_step
+  mass_water_u = mass_water_u_last + prcp_u * length_step * density_water * (1 - exp(-lai_u * clumping))
+  massMax_water_u = 0.1 * lai_u
+
+  mass_water_u = clamp(mass_water_u, 0, massMax_water_u)
+
+  massStep_water_u = mass_water_u - mass_water_u_last
+  massStep_water_u = max(0.0, massStep_water_u)
+
+  percent_water_u = mass_water_u / massMax_water_u
+  percent_water_u = min(1, percent_water_u)
+
+  # ground
+  prcp_g = prcp_u - massStep_water_u / density_water / length_step
+
+  return mass_water_o, mass_water_u, percent_water_o, percent_water_u, prcp_g
+end
+
+
+function rainfall_stage2_jl(evapo_water_o::Float64, evapo_water_u::Float64,
+  mass_water_o::Ref{Float64}, mass_water_u::Ref{Float64})
+
+  length_step = kstep  # 6min or 360s per step
+  mass_water_o[] = max(mass_water_o[] - evapo_water_o * length_step, 0.0)
+  mass_water_u[] = max(mass_water_u[] - evapo_water_u * length_step, 0.0)
+end
diff --git a/src/beps_inter_prg.jl b/src/beps_inter_prg.jl
index 75d5c74..f69e502 100644
--- a/src/beps_inter_prg.jl
+++ b/src/beps_inter_prg.jl
@@ -184,7 +184,6 @@ function inter_prg_jl(
   end
 
   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)
@@ -204,10 +203,8 @@ function inter_prg_jl(
       Ref(var.alpha_v_sw, kkk), Ref(var.alpha_n_sw, kkk))
 
     # /*****  Rain fall stage 1 by X. Luo  *****/
-    rainfall_stage1(Ta, precip, var.Wcl_o[kkk-1], var.Wcl_u[kkk-1],
-      lai_o, lai_u, clumping,
-      Ref(var.Wcl_o, kkk), Ref(var.Wcl_u, kkk), Ref(var.Xcl_o, kkk),
-      Ref(var.Xcl_u, kkk), Ref(var.r_rain_g, kkk))
+    var.Wcl_o[kkk], var.Wcl_u[kkk], var.Xcl_o[kkk], var.Xcl_u[kkk], var.r_rain_g[kkk] = 
+      rainfall_stage1_jl(Ta, precip, var.Wcl_o[kkk-1], var.Wcl_u[kkk-1], lai_o, lai_u, clumping)
 
     # Old version
     # if(thetam[0][kkk-1]<soilp.theta_vwp[1]*0.5) var.alpha_g = alpha_dry;
@@ -330,7 +327,7 @@ function inter_prg_jl(
       var.Xcl_o[kkk], var.Xcl_u[kkk], var.Xcs_o[kkk], var.Xcs_u[kkk])
 
     # /*****  Rainfall stage 2 by X. Luo  *****/
-    rainfall_stage2(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))
 
     # /*****  Snow pack stage 2 by X. Luo  *****/
     snowpack_stage2(var.EiS_o[kkk], var.EiS_u[kkk], Ref(var.Wcs_o, kkk), Ref(var.Wcs_u, kkk))
diff --git a/src/clang/BEPS_c.jl b/src/clang/BEPS_c.jl
index 7205d6c..c315867 100644
--- a/src/clang/BEPS_c.jl
+++ b/src/clang/BEPS_c.jl
@@ -92,25 +92,32 @@ function lai2(clumping, CosZs, stem_o, stem_u, lai_o, lai_u, LAI::Leaf, PAI::Lea
     clumping, CosZs, stem_o, stem_u, lai_o, lai_u, Ref(LAI), Ref(PAI))
 end
 
-function rainfall_stage1(temp_air, precipitation, mass_water_o_last, mass_water_u_last, lai_o, lai_u, clumping, mass_water_o, mass_water_u, percent_water_o, percent_water_u, precipitation_g)
-  ccall((:rainfall_stage1, libbeps), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), temp_air, precipitation, mass_water_o_last, mass_water_u_last, lai_o, lai_u, clumping, mass_water_o, mass_water_u, percent_water_o, percent_water_u, precipitation_g)
+function meteo_pack(temp, rh, meteo_pack_output)
+  ccall((:meteo_pack, libbeps), Cvoid, (Cdouble, Cdouble, Ptr{Cdouble}), temp, rh, meteo_pack_output)
 end
 
-function rainfall_stage2(evapo_water_o, evapo_water_u, mass_water_o, mass_water_u)
-  ccall((:rainfall_stage2, libbeps), Cvoid, (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), evapo_water_o, evapo_water_u, mass_water_o, mass_water_u)
-end
+function rainfall_stage1(temp_air, prcp, mass_water_o_last, mass_water_u_last, lai_o, lai_u, clumping)
+  mass_water_o = init_dbl()
+  mass_water_u = init_dbl()
+  percent_water_o = init_dbl()
+  percent_water_u = init_dbl()
+  prcp_g = init_dbl()
+
+  ccall((:rainfall_stage1, libbeps), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), temp_air, prcp, mass_water_o_last, mass_water_u_last, lai_o, lai_u, clumping, mass_water_o, mass_water_u, percent_water_o, percent_water_u, prcp_g)
 
-# no prototype is found for this function at beps.h:181:6, please use with caution
-function rainfall_stage3()
-  ccall((:rainfall_stage3, libbeps), Cvoid, ())
+  mass_water_o[], mass_water_u[], percent_water_o[], percent_water_u[], prcp_g[]
 end
 
-function meteo_pack(temp, rh, meteo_pack_output)
-  ccall((:meteo_pack, libbeps), Cvoid, (Cdouble, Cdouble, Ptr{Cdouble}), temp, rh, meteo_pack_output)
+function rainfall_stage2(evapo_water_o, evapo_water_u, mass_water_o::Ref, mass_water_u::Ref)
+  ccall((:rainfall_stage2, libbeps), Cvoid, (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), evapo_water_o, evapo_water_u, mass_water_o, mass_water_u)
 end
 
+# # no prototype is found for this function at beps.h:181:6, please use with caution
+# function rainfall_stage3()
+#   ccall((:rainfall_stage3, libbeps), Cvoid, ())
+# end
 
-function snowpack_stage1(temp_air, precipitation,
+function snowpack_stage1(temp_air, prcp,
   mass_snow_o_last, mass_snow_u_last, mass_snow_g_last,
   mass_snow_o::TypeRef, mass_snow_u::TypeRef, mass_snow_g::TypeRef,
   lai_o::T, lai_u::T, clumping::T,
@@ -118,7 +125,7 @@ function snowpack_stage1(temp_air, precipitation,
   density_snow::TypeRef, depth_snow::TypeRef, albedo_v_snow::TypeRef, albedo_n_snow::TypeRef) where {T<:Real}
 
   ccall((:snowpack_stage1, libbeps), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
-    temp_air, precipitation,
+    temp_air, prcp,
     mass_snow_o_last, mass_snow_u_last, mass_snow_g_last,
     mass_snow_o, mass_snow_u, mass_snow_g, # by reference
     lai_o, lai_u, clumping,
@@ -154,7 +161,7 @@ export inter_prg_c,
   # latent_heat!, 
   # Leaf_Temperatures,
   # evaporation_canopy, 
-  rainfall_stage1, rainfall_stage2, rainfall_stage3, snowpack_stage1, snowpack_stage2, snowpack_stage3
+  rainfall_stage1, rainfall_stage2, snowpack_stage1, snowpack_stage2, snowpack_stage3
 
 
 end # module