diff --git a/src/+equations/calc_rssrbs.m b/src/+equations/calc_rssrbs.m
index 6be64a5d..8d4ed93d 100644
--- a/src/+equations/calc_rssrbs.m
+++ b/src/+equations/calc_rssrbs.m
@@ -1,5 +1,13 @@
-function [rss, rbs] = calc_rssrbs(SMC, LAI, rbs)
-    global SaturatedMC ResidualMC fieldMC
+function [rss, rbs] = calc_rssrbs(SMC, LAI, rbs, ResidualMC, fieldMC)
+    %{
+    This function is to calculate surface resistance, following C. van der Tol
+    et al. 2014 (BG). The emperical parameters "aa, 4.1, and 4.3" are
+    user-defined. And these parameters can be calibrated during winter season
+    when LAI is low.
+    %}
     aa = 3.8;
-    rss = exp((aa + 4.1) - aa * (SMC - ResidualMC(1)) / (fieldMC(1) - ResidualMC(1)));
-    rbs            = rbs * LAI / 4.3;
+    bb = 4.1;
+    cc = 4.3;
+
+    rss = exp((aa + bb) - aa * (SMC - ResidualMC(1)) / (fieldMC(1) - ResidualMC(1)));
+    rbs            = rbs * LAI / cc;
diff --git a/src/+io/select_input.m b/src/+io/select_input.m
index b5ac6713..967cebbf 100644
--- a/src/+io/select_input.m
+++ b/src/+io/select_input.m
@@ -1,4 +1,4 @@
-function [SoilProperties, leafbio, canopy, meteo, angles, SpaceTimeInfo] = select_input(ScopeParameters, digitsVector, canopy, options, SpaceTimeInfo, SoilProperties)
+function [SoilProperties, leafbio, canopy, meteo, angles, SpaceTimeInfo] = select_input(ScopeParameters, digitsVector, canopy, options, SoilProperties, SpaceTimeInfo)
     global Theta_LL theta_s0
     SoilProperties.spectrum      = ScopeParameters.spectrum(digitsVector(16));
     SoilProperties.rss           = ScopeParameters.rss(digitsVector(17));
@@ -82,7 +82,7 @@
         SoilProperties.GAM  = equations.Soil_Inertia0(SoilProperties.cs, SoilProperties.rhos, SoilProperties.lambdas);
     end
     if options.calc_rss_rbs
-        [SoilProperties.rss, SoilProperties.rbs] = equations.calc_rssrbs(SoilProperties.SMC, canopy.LAI, SoilProperties.rbs);
+        [SoilProperties.rss, SoilProperties.rbs] = equations.calc_rssrbs(SoilProperties.SMC, canopy.LAI, SoilProperties.rbs, SoilProperties.ResidualMC, SoilProperties.fieldMC);
     end
 
     if leafbio.Type
diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m
index bbdefc04..d0919576 100644
--- a/src/STEMMUS_SCOPE.m
+++ b/src/STEMMUS_SCOPE.m
@@ -212,7 +212,7 @@
 ScopeParametersNames = fieldnames(ScopeParameters);
 if options.simulation == 1
     vi = ones(length(ScopeParametersNames), 1);
-    [soil, leafbio, canopy, meteo, angles, xyt]  = io.select_input(ScopeParameters, vi, canopy, options);
+    [soil, leafbio, canopy, meteo, angles, xyt]  = io.select_input(ScopeParameters, vi, canopy, options, SoilProperties);
     [ScopeParameters, xyt, canopy]  = io.loadTimeSeries(ScopeParameters, leafbio, soil, canopy, meteo, constants, F, xyt, path_input, options);
 else
     soil = struct;
@@ -512,7 +512,7 @@
         if options.simulation == 0
             vi(vmax == telmax) = k;
         end
-        [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, vi, canopy, options, xyt, soil);
+        [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, vi, canopy, options, soil, xyt);
         if options.simulation ~= 1
             fprintf('simulation %i ', k);
             fprintf('of %i \n', telmax);