Potential bug in mo_mpr_soilmoist
I think I found a bug in this loop calculating the soil properties:
do i = 1, size(is_present)
if (is_present(i) .lt. 1) cycle
horizon : do j = 1, nHorizons(i)
! calculating vertical hydraulic conductivity
call hydro_cond(Ks_tmp, param(10 : 13), sand(i, j), clay(i, j))
Ks(i, j, :) = Ks_tmp
! calculating other soil hydraulic properties
! tillage horizons
if (j .le. nTillHorizons(i)) then
! LC class
do L = 1, max_LCover
select case (L)
case(1) ! forest
pOM = tmp_orgMatterContent_forest
case(2) ! impervious
pOM = tmp_orgMatterContent_impervious !param(2)
case(3) ! permeable
pOM = tmp_orgMatterContent_pervious
case default
stop 'Error mpr_sm: pOM used uninitialized.'
end select
pM = 100.0_dp - pOM
! bulk density acording to Rawl's (1982) paper
Db(i, j, L) = 100.0_dp / ((pOM / BulkDens_OrgMatter) + (pM / DbM(i, j)))
! Effect of organic matter content
! This is taken into account in a simplified form by using
! the ratio of(Bd / BdOM)
Ks_tmp = Ks_tmp * (DbM(i, j) / Db(i, j, L)) !!!!!!!!! <-- this is the critical line, Ks_tmp gets updated
Ks(i, j, L) = Ks_tmp
! estimated SMs_till & van Genuchten's shape parameter (n)
call Genuchten(thetaS_till(i, j, L), Genu_Mual_n, Genu_Mual_alpha, &
param(4 : 9), sand(i, j), clay(i, j), Db(i, j, L))
! estimating field capacity
call field_cap(thetaFC_till(i, j, L), Ks_tmp, thetaS_till(i, j, L), Genu_Mual_n)
! estimating permanent wilting point
call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS_till(i, j, L), thetaPW_till(i, j, L))
end do
! deeper layers
else
! estimate SMs & van Genuchten's shape parameter (n)
call Genuchten(thetaS(i, j - tmp_minSoilHorizon), Genu_Mual_n, Genu_Mual_alpha, &
param(4 : 9), sand(i, j), clay(i, j), DbM(i, j))
! estimate field capacity
call field_cap(thetaFC(i, j - tmp_minSoilHorizon), &
Ks_tmp, thetaS(i, j - tmp_minSoilHorizon), Genu_Mual_n)
! estimate permanent wilting point
call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS(i, j - tmp_minSoilHorizon), &
thetaPW(i, j - tmp_minSoilHorizon))
end if
end do horizon
end do
The problem is with the line Ks_tmp = Ks_tmp * (DbM(i, j) / Db(i, j, L))
In each iteration over the the land_covers, Ks_tmp
gets updated. Effectively, the Ks_tmp
for pervious land cover type (3) is Ks_tmp * DbMineral / Dbforest * DbMineral / Dbimpervious * DbMineral / Dbpervious
.