Skip to content
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

Updates in Porous Media and Runoff Properties. #2

Merged
merged 2 commits into from
Dec 17, 2024
Merged

Updates in Porous Media and Runoff Properties. #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8ed0374
Modification of field capacity calculation and parallelization in run…
anaioliveira Dec 16, 2024
141008f
Update Porous Media
anaioliveira Dec 17, 2024
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
128 changes: 82 additions & 46 deletions Software/MOHIDLand/ModulePorousMedia.F90
View file
Open in desktop
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