diff --git a/src/global.f90 b/src/global.f90 index abb27f46..07ef2682 100644 --- a/src/global.f90 +++ b/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 diff --git a/src/simul.f90 b/src/simul.f90 index f200c362..a18ad750 100644 --- a/src/simul.f90 +++ b/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,10 +1246,10 @@ 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) @@ -1252,6 +1257,7 @@ subroutine calculate_transpiration(Tpot, Coeffb0Salt, Coeffb1Salt, Coeffb2Salt) 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