Skip to content

Commit

Permalink
Merge pull request #2 from anaioliveira/dev_ana
Browse files Browse the repository at this point in the history
Updates in Porous Media and Runoff Properties.
  • Loading branch information
anaioliveira authored Dec 17, 2024
2 parents 05975aa + 141008f commit 74d2f38
Show file tree
Hide file tree
Showing 2 changed files with 93 additions and 56 deletions.
128 changes: 82 additions & 46 deletions Software/MOHIDLand/ModulePorousMedia.F90
Original file line number Diff line number Diff line change
Expand Up @@ -830,7 +830,7 @@ subroutine ConstructPorousMedia(ModelName, &
endif

!vegetation model growth needs field capacity computation
if (Me%SoilOpt%ComputeSoilField) then
if (Me%SoilOpt%ComputeSoilField .or. (Me%ObjDrainageNetwork /= 0 .and. Me%SoilOpt%CalcDrainageNetworkFlux)) then !Modified by Ana Oliveira - field capacity is needed for PM-DN fluxes
call ComputeSoilFieldCapacity
endif

Expand Down Expand Up @@ -1150,7 +1150,7 @@ subroutine ReadDataFile
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'COMPUTE_SOIL_FIELD', &
Default = .false., &
Default = .true., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR080'
Expand Down Expand Up @@ -2312,7 +2312,7 @@ subroutine AllocateVariables
allocate(Me%EfectiveEVTP (ILB:IUB,JLB:JUB))
allocate(Me%ImpermeableFraction (ILB:IUB,JLB:JUB))

if (Me%SoilOpt%ComputeSoilField) then
if (Me%SoilOpt%ComputeSoilField .or. (Me%ObjDrainageNetwork /= 0 .and. Me%SoilOpt%CalcDrainageNetworkFlux)) then
allocate (Me%ThetaField (ILB:IUB,JLB:JUB,KLB:KUB))
endif

Expand Down Expand Up @@ -4180,7 +4180,7 @@ subroutine StartWithFieldCapacity

!Local-----------------------------------------------------------------
integer :: i,j,k
real :: Hinf
!real :: Hinf
integer :: WTCell
real :: inf_border

Expand Down Expand Up @@ -4211,14 +4211,17 @@ subroutine StartWithFieldCapacity
!Sets cell above WT to Field Capacity
do k = WTCell, Me%WorkSize%KUB

if (k == WTCell) then
!Half of the cell distance
Hinf = - Me%ExtVar%DWZ(i,j,k) * 0.5
else
Hinf = - (Me%ExtVar%DZZ(i,j,k-1) - Hinf)
endif
!if (k == WTCell) then
! !Half of the cell distance
! Hinf = - Me%ExtVar%DWZ(i,j,k) * 0.5
!else
! Hinf = - (Me%ExtVar%DZZ(i,j,k-1) - Hinf)
!endif
!
!Me%Theta(i, j, k) = Theta_(Hinf, Me%SoilID(i, j, k))

Me%Theta(i, j, k) = Theta_(Hinf, Me%SoilID(i, j, k))
!Uses new function to estimate field capacity - Ana Oliveira
Me%Theta(i, j, k) = ThetaFieldCapacity_(Me%SoilID(i, j, k))

enddo

Expand All @@ -4239,7 +4242,7 @@ subroutine ComputeSoilFieldCapacity

!Local-----------------------------------------------------------------
integer :: i,j,k
real :: Hinf
!real :: Hinf
!Begin-----------------------------------------------------------------

do j= Me%WorkSize%JLB, Me%WorkSize%JUB
Expand All @@ -4248,14 +4251,17 @@ subroutine ComputeSoilFieldCapacity

do k = Me%ExtVar%KFloor(i, j), Me%WorkSize%KUB

if (k == Me%ExtVar%KFloor(i, j)) then
!Half of the cell distance
Hinf = - Me%ExtVar%DWZ(i,j,k) * 0.5
else
Hinf = - (Me%ExtVar%DZZ(i,j,k-1) - Hinf)
endif
!if (k == Me%ExtVar%KFloor(i, j)) then
! !Half of the cell distance
! Hinf = - Me%ExtVar%DWZ(i,j,k) * 0.5
!else
! Hinf = - (Me%ExtVar%DZZ(i,j,k-1) - Hinf)
!endif
!
!Me%ThetaField(i, j, k) = Theta_(Hinf, Me%SoilID(i, j, k))

Me%ThetaField(i, j, k) = Theta_(Hinf, Me%SoilID(i, j, k))
!Uses new function to estimate field capacity - Ana Oliveira
Me%ThetaField(i, j, k) = ThetaFieldCapacity_(Me%SoilID(i, j, k))

enddo

Expand Down Expand Up @@ -6590,6 +6596,8 @@ subroutine VariableSaturatedFlow(InfiltrationColumn)
!remove/add fluxes on all layers that are influenced by level difference
elseif (Me%SoilOpt%DNLink == GWFlowToChanByLayer_) then
call ExchangeWithDrainageNet_3
!elseif (Me%SoilOpt%DNLink == GWFlowToChanByHeadPressure_) then
! call ExchangeWithDrainageNet_4
endif
endif

Expand Down Expand Up @@ -8304,7 +8312,7 @@ subroutine SoilParameters

!Close or below residual value
else if (Me%Theta(i,j,k) < Me%RC%ThetaR(i,j,k) + Me%CV%LimitThetaLo) then

!This creates mass...
Me%RC%ThetaF (i, j, k) = Me%CV%LimitThetaLo / (Me%SoilTypes(Me%SoilID(I,J,K))%ThetaS &
- Me%SoilTypes(Me%SoilID(I,J,K))%ThetaR)
Expand Down Expand Up @@ -8918,11 +8926,11 @@ subroutine CalculateUGWaterLevel
!Local-----------------------------------------------------------------
integer :: i, j, k, CHUNK
real :: CellBottomLevel, DZInCell
real :: FieldHead, FieldTheta
!real :: FieldHead, FieldTheta

CHUNK = ChunkJ !CHUNK_J(Me%WorkSize%JLB, Me%WorkSize%JUB)

!$OMP PARALLEL PRIVATE(I,J,K, CellBottomLevel, DZInCell, FieldHead, FieldTheta)
!$OMP PARALLEL PRIVATE(I,J,K, CellBottomLevel, DZInCell)
!$OMP DO SCHEDULE(DYNAMIC, CHUNK)
do J = Me%WorkSize%JLB, Me%WorkSize%JUB
do I = Me%WorkSize%ILB, Me%WorkSize%IUB
Expand All @@ -8942,13 +8950,13 @@ subroutine CalculateUGWaterLevel
k = min(k, Me%WorkSize%KUB)

CellBottomLevel = - Me%ExtVar%SZZ(i,j,k-1)
FieldHead = - 0.5 * Me%ExtVar%DWZ(i, j, k)
FieldTheta = Theta_(FieldHead, Me%SoilID(i,j,k))
!FieldHead = - 0.5 * Me%ExtVar%DWZ(i, j, k)
!FieldTheta = Theta_(FieldHead, Me%SoilID(i,j,k))

if (Me%Theta(i,j,k) < FieldTheta) then
if (Me%Theta(i,j,k) < Me%ThetaField(i, j, k)) then
DZInCell = 0.0
else
DZInCell = LinearInterpolation(FieldTheta, 0.0, Me%RC%ThetaS(i,j,k), &
DZInCell = LinearInterpolation(Me%ThetaField(i, j, k), 0.0, Me%RC%ThetaS(i,j,k), &
Me%ExtVar%DWZ(i,j,k), Me%Theta(i,j,k))
endif

Expand Down Expand Up @@ -9038,6 +9046,31 @@ real function Theta_ (head, SoilID)

end function Theta_

!----------------------------------------------------------------------------

! Added by Ana Oliveira to correct the calculation of field capacity
real function ThetaFieldCapacity_ (SoilID)

!Arguments-------------------------------------------------------------------
integer, intent(IN) :: SoilID

!----------------------------------------------------------------------
real :: DrainageFlux
!----------------------------------------------------------------------

! Formula from Twarakavi, N. K. C, M. Sakai, and J. �im?nek, An objective analysis
! of the dynamic nature of field capacity, Water Resources Research, 45, W10410,
! doi:10.1029/2009WR007944, 9 pp., 2009.

! Considering a drainage flux of 0.01 cm d-1
! m/s = cm/d m s
DrainageFlux = 0.01 * 1E-2 / 86400
ThetaFieldCapacity_ = Me%SoilTypes(SoilID)%ThetaR + &
Me%SoilTypes(SoilID)%Nfit ** (0.6 * LOG10(DrainageFlux/Me%SoilTypes(SoilID)%SatK)) * &
(Me%SoilTypes(SoilID)%ThetaS - Me%SoilTypes(SoilID)%ThetaR)

end function ThetaFieldCapacity_

!----------------------------------------------------------------------------

subroutine ComputeThetaForHead (head, array, mask)
Expand Down Expand Up @@ -9252,7 +9285,7 @@ subroutine ExchangeWithDrainageNet_1
!Local-----------------------------------------------------------------
integer :: i, j, k, STAT_CALL
real :: dH, dX, TotalArea !, LayerArea
real :: MinHead, FieldTheta !, TopHeightInCell
!real :: MinHead, FieldTheta , TopHeightInCell
real :: InfiltrationVolume, ChannelsVolume, WaveHeight
real :: MaxFlow !, VerticalFluxVariation
real, dimension(:, :), pointer :: ChannelsWaterLevel
Expand Down Expand Up @@ -9282,14 +9315,14 @@ subroutine ExchangeWithDrainageNet_1


!!Computation of TotalFlow - always
!$OMP PARALLEL PRIVATE(I,J,K, dH, dX, TotalArea, MinHead, FieldTheta), &
!$OMP PARALLEL PRIVATE(I,J,K, dH, dX, TotalArea), &
!$OMP& PRIVATE(InfiltrationVolume, ChannelsVolume,MaxFlow)
!$OMP DO SCHEDULE(DYNAMIC, CHUNKJ)
do1: do J = Me%WorkSize%JLB, Me%WorkSize%JUB
do2: do I = Me%WorkSize%ILB, Me%WorkSize%IUB

if (Me%ExtVar%RiverPoints(i, j) == OpenPoint) then
if ((ChannelsWaterLevel(i, j) - ChannelsBottomLevel(i, j)) > 0.0) then
if ((ChannelsWaterLevel(i, j) - ChannelsBottomLevel(i, j)) >= 0.0) then

! if (ChannelsBottomLevel(i,j) < Me%ExtVar%BottomTopoG(i, j)) then
! write (*,*)
Expand Down Expand Up @@ -9343,7 +9376,7 @@ subroutine ExchangeWithDrainageNet_1
! * Me%SoilOpt%FCHCondFactor !Me%SoilOpt%HCondFactor
Me%lFlowToChannels(i, j) = (min(dH, WaveHeight) / dX ) * TotalArea * Me%SatK(i, j, k) &
* Me%SoilOpt%FCHCondFactor

!If the channel looses water (infiltration), then set max flux so that volume in channel does not get
!negative.
if (dH < 0) then
Expand Down Expand Up @@ -9377,16 +9410,16 @@ subroutine ExchangeWithDrainageNet_1
!will be half the height
elseif (dH > 0) then
!k = Me%UGCell(i,j)
MinHead = -0.5 * Me%ExtVar%DWZ(i, j, k)
FieldTheta = Theta_(MinHead, Me%SoilID(i,j,k))
!MinHead = -0.5 * Me%ExtVar%DWZ(i, j, k)
!FieldTheta = Theta_(MinHead, Me%SoilID(i,j,k))

!Maxflow was erroneous. It has to use actual Theta and not ThetaS because content may not be
!completely saturated
!m3/s = (-) * m3 / s
if (Me%Theta (i,j,k) <= FieldTheta) then
if (Me%Theta (i,j,k) <= Me%ThetaField(i, j, k)) then
MaxFlow = 0.
else
MaxFlow = (Me%Theta (i,j,k) - FieldTheta) * Me%ExtVar%CellVolume(i,j,k) / Me%CV%CurrentDT
MaxFlow = (Me%Theta (i,j,k) - Me%ThetaField(i, j, k)) * Me%ExtVar%CellVolume(i,j,k) / Me%CV%CurrentDT
!MaxFlow = (Me%Theta (i,j,k) - FieldTheta) * Me%ExtVar%CellVolume(i,j,k) / Me%ExtVar%DT
!MaxFlow = (Me%RC%ThetaS (i,j,k) - FieldTheta) * Me%ExtVar%CellVolume(i,j,k) / Me%ExtVar%DT
endif
Expand Down Expand Up @@ -9560,7 +9593,7 @@ subroutine ExchangeWithDrainageNet_2
integer :: i, j, k, STAT_CALL, CHUNK
real :: dH, dX, TotalArea
real :: TotalHeight, Toplevel, BottomLevel, ChannelColumn
real :: MinHead, FieldTheta, WaveHeight
real :: WaveHeight !FieldTheta, MinHead
real :: InfiltrationVolume, ChannelsVolume
real :: MaxFlow
real, dimension(:, :), pointer :: ChannelsWaterLevel
Expand Down Expand Up @@ -9591,7 +9624,7 @@ subroutine ExchangeWithDrainageNet_2
CHUNK = ChunkJ !CHUNK_J(Me%WorkSize%JLB, Me%WorkSize%JUB)

!$OMP PARALLEL PRIVATE(I,J,K,dH,ChannelColumn,Toplevel,BottomLevel,TotalHeight,TotalArea,dX), &
!$OMP& PRIVATE(ChannelsVolume,InfiltrationVolume,MinHead,FieldTheta,MaxFlow)
!$OMP& PRIVATE(ChannelsVolume,InfiltrationVolume,MaxFlow)
!$OMP DO SCHEDULE(DYNAMIC, CHUNK)
do1: do J = Me%WorkSize%JLB, Me%WorkSize%JUB
do2: do I = Me%WorkSize%ILB, Me%WorkSize%IUB
Expand Down Expand Up @@ -9667,14 +9700,14 @@ subroutine ExchangeWithDrainageNet_2
!will be half the height
elseif (dH > 0) then
!k = Me%UGCell(i,j)
MinHead = -0.5 * Me%ExtVar%DWZ(i, j, k)
FieldTheta = Theta_(MinHead, Me%SoilID(i,j,k))
!MinHead = -0.5 * Me%ExtVar%DWZ(i, j, k)
!FieldTheta = Theta_(MinHead, Me%SoilID(i,j,k))

!m3/s = (-) * m3 / s
if (Me%Theta (i,j,k) <= FieldTheta) then
if (Me%Theta (i,j,k) <= Me%ThetaField(i, j, k)) then
MaxFlow = 0.
else
MaxFlow = (Me%Theta (i,j,k) - FieldTheta) * Me%ExtVar%CellVolume(i,j,k) / Me%CV%CurrentDT
MaxFlow = (Me%Theta (i,j,k) - Me%ThetaField(i, j, k)) * Me%ExtVar%CellVolume(i,j,k) / Me%CV%CurrentDT
endif

if (Me%lFlowToChannels(i,j) > MaxFlow) then
Expand Down Expand Up @@ -9782,7 +9815,7 @@ subroutine ExchangeWithDrainageNet_3
dH = Me%UGWaterLevel2D(i, j) - ChannelsWaterLevel(i, j)
ChannelColumn = ChannelsWaterLevel(i,j) - ChannelsBottomLevel(i, j)
WaveHeight = max(Me%UGWaterLevel2D(i, j), ChannelsWaterLevel(i, j)) - Me%ExtVar%BottomTopoG(i, j)

!Flux only occurs if there is gradient or water
if (((dH > 0.0) .and. ((Me%UGWaterLevel2D(i, j)-Me%ExtVar%BottomTopoG(i, j))>0.0)) .or. &
(dH < 0 .and. ChannelColumn > 0.)) then
Expand Down Expand Up @@ -9818,7 +9851,7 @@ subroutine ExchangeWithDrainageNet_3

!go trough all active flux layers (computed in the routine GetLayersLimits
do3: do K = Me%FlowToChannelsBottomLayer(i,j), Me%FlowToChannelsTopLayer(i,j)

!if layers equal the height is level difference
if (Me%FlowToChannelsTopLayer(i,j) == Me%FlowToChannelsBottomLayer(i,j)) then

Expand Down Expand Up @@ -9887,9 +9920,13 @@ subroutine ExchangeWithDrainageNet_3
!m3/s = m * m * m / s
MaxFlow = - LayerHeight * ChannelsBottomWidth(i, j) * ChannelsNodeLength(i, j) / Me%CV%CurrentDT
else
!limit to residual
!m3/s = m3H20/m3cell * m3cell / s
MaxFlow = (Me%Theta(i,j,k) - Me%RC%ThetaR(i,j,k)) * Me%ExtVar%CellVolume(i,j,k) / Me%CV%CurrentDT
!limit to field theta or zero if theta is higher than field theta
if (Me%Theta(i,j,k) > Me%ThetaField(i, j, k)) then
!m3/s = m3H20/m3cell * m3cell / s
MaxFlow = (Me%Theta(i,j,k) - Me%ThetaField(i, j, k)) * Me%ExtVar%CellVolume(i,j,k) / Me%CV%CurrentDT
else
MaxFlow = 0
endif
endif

if (abs(Me%lFlowToChannelsLayer(i, j, k)) > abs(MaxFlow)) then
Expand Down Expand Up @@ -9956,7 +9993,6 @@ subroutine ExchangeWithDrainageNet_3
end subroutine ExchangeWithDrainageNet_3

!--------------------------------------------------------------------------

subroutine GetLayersLimits_2

!Arguments-------------------------------------------------------------
Expand Down
Loading

0 comments on commit 74d2f38

Please sign in to comment.