From 555cf2e4d23460ed45827f38f79f6f228e8f6c7c Mon Sep 17 00:00:00 2001
From: DanyangYu <Yudanyang123@gmail.com>
Date: Wed, 24 May 2023 15:45:14 +0200
Subject: [PATCH] fix issue #163

---
 src/Root_properties.m | 59 ++++++++++++++++++-------------------------
 src/ebal.m            |  2 +-
 2 files changed, 26 insertions(+), 35 deletions(-)

diff --git a/src/Root_properties.m b/src/Root_properties.m
index 27130beb..d2659866 100644
--- a/src/Root_properties.m
+++ b/src/Root_properties.m
@@ -3,43 +3,27 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% REFERENCES
-function [Rl] = Root_properties(Rl, Ac, rroot, frac, bbx, KT)
-    %%% INPUTS
-    global DeltZ sfactor LAI_msr
-    % BR = 10:1:650; %% [gC /m^2 PFT]
-    % rroot =  0.5*1e-3 ; % 3.3*1e-4 ;%% [0.5-6 *10^-3] [m] root radius
-    %%% OUTPUTS
-    if KT < 2880
-        fr = 0.3 * 3 * exp(-0.15 * LAI_msr(KT)) / (exp(-0.15 * LAI_msr(KT)) + 2 * sfactor);
-        if fr < 0.15
-            fr = 0.15;
-        end
-    else
-        fr = 0.001;
-    end
-    % if KT<3840
-    % fr=-0.091*LAI(KT)+0.40;
-    % else
-    % fr=0.001;
-    % end
-    % elseif KT<3200
-    % fr=0.05;
-    % else
-    % fr=0.02;
-    % end
+%{
+    This function is used to calculate the dynamic growth of root
+    This part can refer to Wang et al. (2021) "Intergrated modelling of photosynthesis and transfer of energy, mass, and momentum in the SPAC system"
+%}
+
+function [Rl] = Root_properties(Rl, Ac, rroot, frac, bbx, KT, DeltZ, sfactor, LAI_msr)
+%     %%% INPUTS
+%     global DeltZ sfactor LAI_msr
+
+    fr = calculateRootfraction(KT);
+
     DeltZ0 = DeltZ' / 100;
     BR = Ac * fr * 1800 * 12 / 1000000;
     root_den = 250 * 1000; %% [gDM / m^3] Root density  Jackson et al., 1997
     R_C = 0.488; %% [gC/gDM] Ratio Carbon-Dry Matter in root   Jackson et al.,  1997
     nn = numel(Rl);
+    
+    % This is used to simulate the root growth
     if (~isnan(Ac)) || (Ac > 0)
         Rl = Rl .* DeltZ0;
         Delta_Rltot = BR / R_C / root_den / (pi * (rroot^2)); %% %% root length index [m root / m^2 PFT]
-        % for i=1:nn
-        %   if Rl(i)>50000
-        %       frac(i)=0;
-        %    end
-        % end
         if ~isnan(frac)
             frac = frac / sum(sum(frac));
         else
@@ -48,11 +32,18 @@
         Delta_Rl = frac .* bbx * Delta_Rltot;
         Rl = Rl + Delta_Rl;
         Rl = Rl ./ DeltZ0;
+    end
+
+end
+
+function fr = calculateRootfraction(KT)
+    % this function is used to calculate the root fraction
+    if KT < 2880 %2880 means the time step when the root stops growing
+        fr = 0.3 * 3 * exp(-0.15 * LAI_msr(KT)) / (exp(-0.15 * LAI_msr(KT)) + 2 * sfactor);
+        if fr < 0.15
+            fr = 0.15;
+        end
     else
+        fr = 0.001;
     end
-    % for i=1:nn
-    %    if Rl(i)>50000
-    %       Rl(i)=50000;
-    %    end
-    % end
 end
diff --git a/src/ebal.m b/src/ebal.m
index 76e07cd6..e1e5cf48 100644
--- a/src/ebal.m
+++ b/src/ebal.m
@@ -461,7 +461,7 @@
     Rntot_PAR       = LAI * (Fc * Rnh_PAR  + equations.meanleaf(canopy, Rnu_PAR, 'angles_and_layers', Ps)); % net PAR leaves
     aPAR_Cab_eta        = LAI * (Fc * (profiles.etah .* Rnh_PAR) + equations.meanleaf(canopy, profiles.etau .* Rnu_PAR, 'angles_and_layers', Ps));
     % ... green ePAR * relative fluorescence emission efficiency
-    %%%%%%%%%%%%%%%%%%% [Delta_Rltot] = Root_properties(Actot,rroot);
+    %%%%%%%%%%%%%%%%%%% [Delta_Rltot] = Root_properties(Rl, Ac, rroot, frac, bbx, KT, DeltZ, sfactor, LAI_msr);
     %%%%%%%%%%%%%%%%%%% Delta_Rl = fc*Delta_Rltot;
     %%%%%%%%%%%%%%%%%%% Rl = Rl + Delta_Rl;
     %%%%%%%%%%%%%%%%%%% Rltot = sum(sum(Rl));