V73: Improved root water extraction profile

src/global.f90

@@ -198,6 +198,10 @@ module ac_global
         !! - relative wetness (RUNOFF)
         !! - evaporation process
         !! - transpiration process *)
+    real(dp) :: SinkMajor 
+        !! maximum possible water extraction (with current water stress)
+    real(dp) :: SinkMinor 
+        !! required extraction considering root distribution (no water stress)
     !! salinity factors
     real(dp), dimension(11) :: Salt
         !! salt content in solution in cells (g/m2)
@@ -15098,6 +15102,24 @@ function GetCompartment_WFactor(i) result(WFactor)
 end function GetCompartment_WFactor
 
 
+function GetCompartment_SinkMajor(i) result(SinkMajor)
+    !! Getter for the "SinkMajor" attribute of the "compartment" global variable.
+    integer(int32), intent(in) :: i
+    real(dp) :: SinkMajor
+
+    SinkMajor = compartment(i)%SinkMajor
+end function GetCompartment_SinkMajor
+
+
+function GetCompartment_SinkMinor(i) result(SinkMinor)
+    !! Getter for the "SinkMinor" attribute of the "compartment" global variable.
+    integer(int32), intent(in) :: i
+    real(dp) :: SinkMinor
+
+    SinkMinor = compartment(i)%SinkMinor
+end function GetCompartment_SinkMinor
+
+
 function GetCompartment_Salt(i1, i2) result(Salt)
     !! Getter for individual elements of "Salt" attribute of the "compartment" global variable.
     integer(int32), intent(in) :: i1
@@ -15199,6 +15221,24 @@ subroutine SetCompartment_WFactor(i, WFactor)
 end subroutine SetCompartment_WFactor
 
 
+subroutine SetCompartment_SinkMajor(i, SinkMajor)
+    !! Setter for the "SinkMajor" attribute of the "compartment" global variable.
+    integer(int32), intent(in) :: i
+    real(dp), intent(in) :: SinkMajor
+
+    compartment(i)%SinkMajor = SinkMajor
+end subroutine SetCompartment_SinkMajor
+
+
+subroutine SetCompartment_SinkMinor(i, SinkMinor)
+    !! Setter for the "SinkMinor" attribute of the "compartment" global variable.
+    integer(int32), intent(in) :: i
+    real(dp), intent(in) :: SinkMinor
+
+    compartment(i)%SinkMinor = SinkMinor
+end subroutine SetCompartment_SinkMinor
+
+
 subroutine SetCompartment_Salt(i1, i2, Salt)
     !! Setter for individual elements of "Salt" attribute of the "compartment global variable.
     integer(int32), intent(in) :: i1

src/simul.f90

@@ -45,6 +45,8 @@ module ac_simul
                      GetCompartment_theta, &
                      GetCompartment_Thickness, &
                      GetCompartment_WFactor, &
+                     GetCompartment_SinkMajor, &
+                     GetCompartment_SinkMinor, &
                      GetCrop, &
                      GetCrop, &
                      GetCrop_aCoeff, &
@@ -294,6 +296,8 @@ module ac_simul
                      SetCompartment_Salt, &
                      SetCompartment_theta, &
                      SetCompartment_WFactor, &
+                     SetCompartment_SinkMajor, &
+                     SetCompartment_SinkMinor, &
                      SetCrop_CCoAdjusted, &
                      SetCrop_CCxAdjusted, &
                      SetCrop_CCxWithered, &
@@ -1166,6 +1170,7 @@ subroutine calculate_transpiration(Tpot, Coeffb0Salt, Coeffb1Salt, Coeffb2Salt)
     real(dp) :: TpotMAX, RedFact, RedFactECsw
     real(dp) :: Wrel, WrelSalt, pStomatLLAct, crop_pActStom_tmp
     real(dp) :: CompiECe, CompiECsw, CompiECswFC
+    real(dp) :: TheSum, TheRatio, TheResidue
     logical :: SWCtopSoilConsidered_temp
     type(CompartmentIndividual), dimension(max_No_compartments) :: Comp_temp
     type(CompartmentIndividual) :: Compi_temp
@@ -1241,17 +1246,18 @@ subroutine calculate_transpiration(Tpot, Coeffb0Salt, Coeffb1Salt, Coeffb2Salt)
                                               GetRootZoneWC_Actual(), &
                                               real(GetCrop_AnaeroPoint(), kind=dp), &
                                               GetRootingDepth(), RedFact)
-            TpotMAX = RedFact * TpotMax
+            TpotMAX = RedFact * TpotMAX
         end if
 
-        ! 2. extraction of TpotMax out of the compartments
+        ! 2. extraction of TpotMAX out of the compartments
         ! 2.a initial settings
         Comp_temp = GetCompartment()
         call calculate_rootfraction_compartment(GetRootingDepth(), Comp_temp)
         call calculate_sink_values(TpotMAX, GetRootingDepth(), Comp_temp, GetCrop())
         call SetCompartment(Comp_temp)
         compi = 0
         pre_layer = 0
+        TheSum = 0._dp
         loop: do
             compi = compi + 1
             layeri = GetCompartment_Layer(compi)
@@ -1307,22 +1313,65 @@ subroutine calculate_transpiration(Tpot, Coeffb0Salt, Coeffb1Salt, Coeffb2Salt)
                 call Correction_Anaeroby(Compi_temp, alfa)
                 call SetCompartment_i(compi, Compi_temp)
             end if
-            ! 2.c extract water
-            sinkMM = 1000._dp * (alfa * GetCompartment_WFactor(compi) * &
-                    GetCompartment_Smax(compi)) * GetCompartment_Thickness(compi)
-            WtoExtract = TpotMAX-GetTact()
-            if (WtoExtract < sinkMM) then
-                sinkMM = WtoExtract
+            ! 2.c determine maximum extration amount (mm) in current root zone
+            ! maximum possible (with current water stress) 
+            call SetCompartment_SinkMajor(compi, 1000 * (alfa * GetCompartment_WFactor(compi) * &
+                GetCompartment_Smax(compi)) * GetCompartment_Thickness(compi))
+            ! Sum of maximum possible (without any water stress)
+            if (GetCompartment_WFactor(compi) > 0.00001) then
+                TheSum = TheSum + 1000 * (GetCompartment_WFactor(compi) * GetCompartment_Smax(compi)) * &
+                         GetCompartment_Thickness(compi)
             end if
-            call SetCompartment_theta(compi, GetCompartment_theta(compi) &
-                 - sinkMM/(1000._dp*GetCompartment_Thickness(compi)* &
-                 (1._dp - GetSoilLayer_GravelVol(layeri)/100._dp)))
-            WtoExtract = WtoExtract - sinkMM
-            call SetTact(GetTact() + sinkMM)
-            if ((WtoExtract < epsilon(1._dp) .or. &
-                                    (compi == GetNrCompartments()))) exit loop
+
+            if (compi == GetNrCompartments()) exit loop
         end do loop
 
+        ! 2.d Determine required extraction considering root distribution
+        TheRatio = TpotMAX / TheSum
+
+        do compi = 1, GetNrcompartments()
+            if (GetCompartment_WFactor(compi) > 1.0E-5) then
+                call SetCompartment_SinkMinor(compi, TheRatio * 1000.0 * &
+                     GetCompartment_WFactor(compi) * GetCompartment_Smax(compi) * &
+                     GetCompartment_Thickness(compi))
+            else
+                call SetCompartment_SinkMinor(compi, 0._dp)
+            end if
+        end do
+
+        ! 2.e Extract water
+        TheResidue = 0.0
+        WtoExtract = TpotMAX
+
+        do compi = 1, GetNrcompartments()
+            if ((GetCompartment_SinkMinor(compi) > 1.0E-5) .and. (WtoExtract > 1.0E-5)) then
+                if (GetCompartment_SinkMinor(compi) < GetCompartment_SinkMajor(compi)) then
+                    if ((GetCompartment_SinkMinor(compi) + TheResidue) > GetCompartment_SinkMajor(compi)) then
+                        TheResidue = TheResidue - (GetCompartment_SinkMajor(compi) - GetCompartment_SinkMinor(compi))
+                        sinkMM = GetCompartment_SinkMajor(compi)
+                    else
+                        sinkMM = GetCompartment_SinkMinor(compi) + TheResidue
+                        TheResidue = 0.0
+                    end if
+                else
+                    sinkMM = GetCompartment_SinkMajor(compi)
+                    TheResidue = TheResidue + (GetCompartment_SinkMinor(compi) - GetCompartment_SinkMajor(compi))
+                end if
+
+                ! Extract water
+                WtoExtract = TpotMAX - GetTact()
+                if (WtoExtract < sinkMM) then
+                    sinkMM = WtoExtract
+                end if
+
+                call SetCompartment_theta(compi, GetCompartment_theta(compi) - &
+                     sinkMM / (1000._dp * GetCompartment_Thickness(compi) * &
+                     (1.0 - GetSoilLayer_GravelVol(layeri) / 100._dp)))
+
+                WtoExtract = WtoExtract - sinkMM
+                call SetTact(GetTact() + sinkMM)
+            end if
+        end do
 
         ! 3. add net irrigation water requirement
         if (GetIrriMode() == IrriMode_Inet) then