Skip to content

V73: Improved root zone expansion #374

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mbechtold wants to merge 5 commits into
base: main
Choose a base branch
from
Open

V73: Improved root zone expansion #374

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1ba05fe
ENH: improved root zone expansion
Jan 23, 2026
f722013
BUG: div by integer
mbechtold Jan 26, 2026
da72708
BUG: For consistency in compartment thickness between Fortran and Pascal
mbechtold Jan 26, 2026
8f21d07
STY: indentation
mbechtold Jan 28, 2026
fddb58d
STY: small value root zone discretization consistency
mbechtold Jan 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 11 additions & 8 deletions src/global.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6606,13 +6606,16 @@ end subroutine LoadInitialConditions
subroutine AdjustSizeCompartments(CropZx)
real(dp), intent(in) :: CropZx

real(dp) :: CropZx_eff
integer(int32) :: i , compi
real(dp) :: TotDepthC, fAdd
integer(int8) :: PrevNrComp
real(dp), dimension(max_No_Compartments) :: PrevThickComp, &
PrevVolPrComp, &
PrevECdSComp

! esnures consistency for adjusted compartment sizes between Pascal and Fortran
CropZx_eff = CropZx + ac_zero_threshold
Copy link
Collaborator

@gdelannoy gdelannoy Jan 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(1) not quite understood why a small value needs to be consistently added in one direction (it is also not really serving as a zero threshold in the same sense as used earlier in the code) (2) typo "ensures"

Copy link
Collaborator

@lbusschaert lbusschaert Jan 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had to do a similar trick in LDT to get the right number of compartments. I think it has basically the same result but instead of having a new CropZx_eff, I adjusted one if statement.
See L185 https://github.com/KUL-RSDA/LISF/blob/working/ac72_GMD/ldt/params/AquaCrop/define_AC_compartments.F90
I am not sure why this happens (small rounding errors, differences in compilation, optimization?), but this line was indeed dangerous
if (TotDepthC < CropZx) then
and is solved by adding a small value

Copy link
Collaborator Author

@mbechtold mbechtold Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, we had the problem earlier. It's because the user defines round values for the max root zone depth like 1.5. Such round values are directly at the edge when the compartment thickness calculation decides to add or remove a 0.05 increment to the total depth of the simulated soil depth. We struggled with that difference in how the Delphi and Pascal code eventually discretizes. Adding a very small value solved the problem. The epsilon was not enough, so I used the ac_zero_threshold which worked.
@gdelannoy It is correct that the meaning of ac_zero_threshold is a bit different. I could also hardcode the small number instead of using this parameter. Should I?

! 1. Save intial soil water profile (required when initial soil
! water profile is NOT reset at start simulation - see 7.)
PrevNrComp = int(GetNrCompartments(), kind=int8)
Expand All @@ -6631,24 +6634,24 @@ subroutine AdjustSizeCompartments(CropZx)
if (GetNrCompartments() < 12) then
loop: do
call SetNrCompartments(GetNrCompartments() + 1)
if ((CropZx - TotDepthC) > GetSimulParam_CompDefThick()) then
if ((CropZx_eff - TotDepthC) > GetSimulParam_CompDefThick()) then
call SetCompartment_Thickness(GetNrCompartments(), &
GetSimulParam_CompDefThick())
else
call SetCompartment_Thickness(GetNrCompartments(), &
CropZx - TotDepthC)
CropZx_eff - TotDepthC)
end if
TotDepthC = TotDepthC &
+ GetCompartment_Thickness(GetNrCompartments())
if ((GetNrCompartments() == max_No_compartments) &
.or. ((TotDepthC + 0.00001) >= CropZx)) exit loop
.or. ((TotDepthC + 0.00001) >= CropZx_eff)) exit loop
end do loop
end if

! 4. Adjust size of compartments (if total depth of compartments < rooting depth)
if ((TotDepthC + 0.00001_dp) < CropZx) then
if ((TotDepthC + 0.00001_dp) < CropZx_eff) then
call SetNrCompartments(12)
fAdd = (CropZx/0.1_dp - 12._dp)/78._dp
fAdd = (CropZx_eff/0.1_dp - 12._dp)/78._dp
do i = 1, 12
call SetCompartment_Thickness(i, 0.1 * (1._dp + i*fAdd))
call SetCompartment_Thickness(i, 0.05 &
Expand All @@ -6659,16 +6662,16 @@ subroutine AdjustSizeCompartments(CropZx)
do i = 1, GetNrCompartments()
TotDepthC = TotDepthC + GetCompartment_Thickness(i)
end do
if (TotDepthC < CropZx) then
if (TotDepthC < CropZx_eff) then
loop2: do
call SetCompartment_Thickness(12, &
GetCompartment_Thickness(12) &
+ 0.05_dp)
TotDepthC = TotDepthC + 0.05_dp
if (TotDepthC >= CropZx) exit loop2
if (TotDepthC >= CropZx_eff) exit loop2
end do loop2
else
do while ((TotDepthC - 0.04999999_dp) >= CropZx)
do while ((TotDepthC - 0.04999999_dp) >= CropZx_eff)
call SetCompartment_Thickness(12, &
GetCompartment_Thickness(12) &
- 0.05_dp)
Expand Down
146 changes: 134 additions & 12 deletions src/run.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,7 @@ module ac_run
GetSaltInfiltr, &
GetSimulation_FromDayNr, &
GetSimulation_IrriECw, &
GetSimulParam_RootPercentZmin, &
GetSimulation_SalinityConsidered, &
GetSimulation_Storage_Btotal, &
GetSimulation_SumGDD, &
Expand Down Expand Up @@ -599,6 +600,8 @@ module ac_run
real(dp) :: BprevSum, YprevSum, SumGDDcuts, HItimesBEF
real(dp) :: ScorAT1, ScorAT2, HItimesAT1, HItimesAT2, HItimesAT
real(dp) :: alfaHI, alfaHIAdj
real(dp) :: tDaysZmin ! time to reach Zmin (Days)
real(dp) :: tGDDZmin ! time to reach Zmin (GDDays)

!! DelayedGermination
integer(int32) :: NextSimFromDayNr !! the Simulation.FromDayNr for next run if delayed germination and KeepSWC
Expand Down Expand Up @@ -5328,6 +5331,11 @@ subroutine InitializeSimulationRunPart2()
GetCrop_RootShape(),&
GetCrop_ModeCycle()))
end if
! 16.3 Time to reach Zmin
call CalculateTimeToReachZmin(GetCrop_RootMin(),GetCrop_RootMax(),GetCrop_RootShape(), &
GetCrop_DaysToGermination(),GetCrop_DaysToMaxRooting(), &
GetCrop_GDDaysToGermination(),GetCrop_GDDaysToMaxRooting(), &
GetCrop_ModeCycle(),tDaysZmin,tGDDZmin)

! 17. Multiple cuttings
call SetNrCut(0)
Expand Down Expand Up @@ -6509,6 +6517,116 @@ subroutine GetPotValSF(DAP, SumGDDAdjCC, PotValSF)
PotValSF = 100._dp * (1._dp/GetCCxCropWeedsNoSFstress()) * PotValSF
end subroutine GetPotValSF


subroutine CalculateTimeToReachZmin(ZMin, ZMax, RootShape, L0, LZmax, GGDL0, GDDLZmax, &
TheModeCycle, tDaysZmin, tGDDZmin)
real(dp), intent(in) :: ZMin
real(dp), intent(in) :: ZMax
integer(int8), intent(in) :: RootShape
integer(int32), intent(in) :: L0
integer(int32), intent(in) :: LZmax
integer(int32), intent(in) :: GGDL0
integer(int32), intent(in) :: GDDLZmax
integer(intEnum), intent(in) :: TheModeCycle
real(dp), intent(inout) :: tDaysZmin
real(dp), intent(inout) :: tGDDZmin

real(dp) :: Zini, Zfunction

tDaysZmin = real(L0, kind=dp)
tGDDZmin = real(GGDL0, kind=dp)

if ((GetSimulParam_RootPercentZmin() < 100) .and. (ZMin < ZMax)) then
Copy link
Collaborator

@gdelannoy gdelannoy Jan 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

100._dp? Not sure anymore if that was needed for >, < operations.

Copy link
Collaborator

@lbusschaert lbusschaert Jan 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's best if we put dp I think (and also consistent with the rest of the code).

Copy link
Collaborator Author

@mbechtold mbechtold Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would say so too. No harm.

Zini = ZMin * (real(GetSimulParam_RootPercentZmin(), kind=dp)/100._dp)

Zfunction = exp( (real(RootShape, kind=dp)/10._dp) * &
log( (ZMin - Zini) / (ZMax - Zini) ) )

select case (TheModeCycle)
case (modeCycle_GDDays)
if (GDDLZmax > (GGDL0/2._dp)) then
tGDDZmin = Zfunction * real((GDDLZmax - (GGDL0/2._dp)), kind=dp) + &
real((GGDL0/2._dp), kind=dp)
end if
case default
if (LZmax > (L0/2._dp)) then
tDaysZmin = Zfunction * real((LZmax - (L0/2._dp)), kind=dp) + &
real((L0/2._dp), kind=dp)
end if
end select
end if
end subroutine CalculateTimeToReachZmin


subroutine CalculateRootingDepth(tDaysZmin, tGDDZmin, ZiPrev, GDDayi, RootingDepth)
real(dp), intent(in) :: tDaysZmin
real(dp), intent(in) :: tGDDZmin
real(dp), intent(in) :: ZiPrev
real(dp), intent(in) :: GDDayi
real(dp), intent(inout) :: RootingDepth

real(dp) :: tCalDayReal, tCalGDD, CalSumGDDPrev
real(dp) :: Zini, Zfunction
integer(int32) :: tCalDay
integer(int32) :: tCalDayNow

! initialize
tCalDayNow = GetDayNri() - GetCrop_Day1() + 1 ! actual time (calendar day)
tCalDay = tCalDayNow
tCalGDD = GetSimulation_SumGDD() ! actual time (sum of GDD)
CalSumGDDPrev = GetSumGDDPrev() ! sum of GDD at previous day

! Adjust time after planting (tCalDay OR tCalGDD) to simulate root zone expansion
Zini = GetCrop_RootMin() * (real(GetSimulParam_RootPercentZmin(), kind=dp)/100._dp)

if ((GetCrop_RootMin() < GetCrop_RootMax()) .and. (ZiPrev > Zini)) then
Copy link
Collaborator

@lbusschaert lbusschaert Jan 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just wondering, is it even possible to have RootMin>RootMax? Shouldn't the code raise an error if it's defined like that in the crop parameters? Now if it does not enter the if statement, nothing happens

Copy link
Collaborator Author

@mbechtold mbechtold Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's probably there to avoid that this part is executed when RootMin == RootMax which then causes errors. True that RootMin should never be higher than RootMax, however, I don't think the current AquaCrop standalone catches non-meaningful parameter combinations elsewhere in the code. The GUI maybe.

! valid for calendar and GDD time
Zfunction = exp( (real(GetCrop_RootShape(), kind=dp)/10._dp) * &
log( (ZiPrev - Zini) / (GetCrop_RootMax() - Zini) ) )

select case (GetCrop_ModeCycle())
case (modeCycle_GDDays)
! growing degree days
if (GetSimulation_SumGDD() > tGDDZmin) then
! root zone expansion can exceed Zmin
! time to simulate root zone expansion is given by Zr of previous day (ZiPrev)
CalSumGDDPrev = Zfunction * &
real((GetCrop_GDDaysToMaxRooting() - (GetCrop_GDDaysToGermination()/2._dp)), kind=dp) + &
real((GetCrop_GDDaysToGermination()/2._dp), kind=dp)

tCalGDD = CalSumGDDPrev + GDDayi ! add GDD of today to get adjusted sum GDD
if (tCalGDD > GetSimulation_SumGDD()) then
tCalGDD = GetSimulation_SumGDD()
end if
end if

case default
! calendar days
if (real(tCalDayNow, kind=dp) > tDaysZmin) then
! root zone expansion can exceed Zmin
! time to simulate root zone expansion is given by Zr of previous day (ZiPrev)
tCalDayReal = Zfunction * &
real((GetCrop_DaysToMaxRooting() - (GetCrop_DaysToGermination()/2._dp)), kind=dp) + &
real((GetCrop_DaysToGermination()/2._dp), kind=dp)

tCalDay = roundc(tCalDayReal, mold=1_int32) + 1 ! integer value + 1 day
if (tCalDay > tCalDayNow) tCalDay = tCalDayNow
end if
end select
end if
! calculate root zone expansion with time after planting given by Zr of previous day (ZiPrev)
RootingDepth = AdjustedRootingDepth( &
GetPlotVarCrop_ActVal(), GetPlotVarCrop_PotVal(), GetTpot(), GetTact(), GetStressLeaf(), &
GetStressSenescence(), tCalDay, GetCrop_DaysToGermination(), GetCrop_DaysToMaxRooting(), GetCrop_DaysToHarvest(), &
GetCrop_GDDaysToGermination(), GetCrop_GDDaysToMaxRooting(), GetCrop_GDDaysToHarvest(), &
CalSumGDDPrev, tCalGDD, &
GetCrop_RootMin(), GetCrop_RootMax(), ZiPrev, GetCrop_RootShape(), &
GetCrop_ModeCycle() )

if (RootingDepth < ZiPrev) RootingDepth = ZiPrev
end subroutine CalculateRootingDepth


!! ===END Subroutines and functions for AdvanceOneTimeStep ===


Expand Down Expand Up @@ -6792,6 +6910,7 @@ subroutine AdvanceOneTimeStep(WPi, HarvestNow)
ScorAT1_temp, ScorAT2_temp, HItimesAT1_temp, &
HItimesAT2_temp, HItimesAT_temp, alfaHI_temp, &
alfaHIAdj_temp, TESTVAL
real(dp) :: RootingDepth_temp
logical :: WaterTableInProfile_temp, NoMoreCrop_temp

! 1. Get ETo
Expand Down Expand Up @@ -6955,18 +7074,21 @@ subroutine AdvanceOneTimeStep(WPi, HarvestNow)
if (((GetDayNri()-GetSimulation_DelayedDays()) >= GetCrop_Day1()) .and. &
((GetDayNri()-GetSimulation_DelayedDays()) <= GetCrop_DayN())) then
! rooting depth at DAP (at Crop.Day1, DAP = 1)
call SetRootingDepth(AdjustedRootingDepth(GetPlotVarCrop_ActVal(), &
GetPlotVarCrop_PotVal(), GetTpot(), GetTact(), &
GetStressLeaf(), GetStressSenescence(), &
(GetDayNri()-GetCrop_Day1()+1), &
GetCrop_DaysToGermination(), &
GetCrop_DaysToMaxRooting(), GetCrop_DaysToHarvest(), &
GetCrop_GDDaysToGermination(), &
GetCrop_GDDaysToMaxRooting(), &
GetCrop_GDDaysToHarvest(), GetSumGDDPrev(), &
(GetSimulation_SumGDD()), GetCrop_RootMin(), &
GetCrop_RootMax(), GetZiprev(), GetCrop_RootShape(), &
GetCrop_ModeCycle()))
call CalculateRootingDepth(tDaysZmin,tGDDZmin,&
GetZiPrev(),GetGDDayi(),RootingDepth_temp)
call SetRootingDepth(RootingDepth_temp)
!call SetRootingDepth(AdjustedRootingDepth(GetPlotVarCrop_ActVal(), &
Copy link
Collaborator

@lbusschaert lbusschaert Jan 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Leave commented? Or remove to clean up?

Copy link
Collaborator Author

@mbechtold mbechtold Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cleaned up

! GetPlotVarCrop_PotVal(), GetTpot(), GetTact(), &
! GetStressLeaf(), GetStressSenescence(), &
! (GetDayNri()-GetCrop_Day1()+1), &
! GetCrop_DaysToGermination(), &
! GetCrop_DaysToMaxRooting(), GetCrop_DaysToHarvest(), &
! GetCrop_GDDaysToGermination(), &
! GetCrop_GDDaysToMaxRooting(), &
! GetCrop_GDDaysToHarvest(), GetSumGDDPrev(), &
! (GetSimulation_SumGDD()), GetCrop_RootMin(), &
! GetCrop_RootMax(), GetZiprev(), GetCrop_RootShape(), &
! GetCrop_ModeCycle()))
call SetZiprev(GetRootingDepth())
! IN CASE rootzone drops below groundwate table
if ((GetZiAqua() >= 0._dp) .and. (GetRootingDepth() > &
Expand Down