Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
12016 lines (8963 sloc) 581 KB
!------------------------------------------------------------------------------
! IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE : Mohid Model
! PROJECT : Mohid Base 1
! MODULE : RunOff
! URL : http://www.mohid.com
! AFFILIATION : IST/MARETEC, Marine Modelling Group
! DATE : Jan 2004
! REVISION : Frank Braunschweig - v4.0
! DESCRIPTION : Module which calculates the Surface RunOff
!
!------------------------------------------------------------------------------
!
!This program is free software; you can redistribute it and/or
!modify it under the terms of the GNU General Public License
!version 2, as published by the Free Software Foundation.
!
!This program is distributed in the hope that it will be useful,
!but WITHOUT ANY WARRANTY; without even the implied warranty of
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!GNU General Public License for more details.
!
!You should have received a copy of the GNU General Public License
!along with this program; if not, write to the Free Software
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!------------------------------------------------------------------------------
Module ModuleRunOff
use ModuleGlobalData
use ModuleTime
use ModuleTimeSerie ,only : StartTimeSerieInput, KillTimeSerie, &
GetTimeSerieInitialData, GetTimeSerieValue, &
StartTimeSerie, WriteTimeSerieLine, &
GetNumberOfTimeSeries, GetTimeSerieLocation, &
TryIgnoreTimeSerie, CorrectsCellsTimeSerie, &
GetTimeSerieName, WriteTimeSerie
use ModuleEnterData
use ModuleHDF5
use ModuleFunctions ,only : TimeToString, SetMatrixValue, ChangeSuffix, &
CHUNK_J, LinearInterpolation, InterpolateValueInTime
use ModuleHorizontalGrid ,only : GetHorizontalGridSize, GetHorizontalGrid, &
UnGetHorizontalGrid, WriteHorizontalGrid, &
GetGridCellArea, GetXYCellZ, &
GetCellZInterceptByLine, &
GetCellZInterceptByPolygon
use ModuleHorizontalMap ,only : GetBoundaries, UngetHorizontalMap
use ModuleGridData ,only : GetGridData, UngetGridData, WriteGridData
use ModuleBasinGeometry ,only : GetBasinPoints, GetRiverPoints, GetCellSlope, &
GetDrainageDirection, TargetPoint, &
UnGetBasin
use ModuleStopWatch ,only : StartWatch, StopWatch
use ModuleFillMatrix ,only : ConstructFillMatrix, ModifyFillMatrix, &
KillFillMatrix
use ModuleDrainageNetwork ,only : GetChannelsWaterLevel, GetChannelsSurfaceWidth, &
GetChannelsBankSlope, GetChannelsNodeLength, &
GetChannelsBottomLevel, UnGetDrainageNetwork, &
GetChannelsID,GetChannelsVolume, &
GetChannelsMaxVolume, GetChannelsActiveState, &
GetChannelsTopArea, GetChannelsVelocity
use ModuleDischarges ,only : Construct_Discharges, GetDischargesNumber, &
GetDischargesGridLocalization, &
GetDischargeWaterFlow, GetDischargesIDName, &
TryIgnoreDischarge, GetDischargeSpatialEmission, &
CorrectsCellsDischarges, Kill_Discharges, &
GetByPassON, GetDischargeFlowDistribuiton, &
UnGetDischarges, SetLocationCellsZ, &
CorrectsBypassCellsDischarges
use ModuleBoxDif, only : StartBoxDif, GetBoxes, GetNumberOfBoxes, UngetBoxDif, &
BoxDif, KillBoxDif
use ModuleDrawing
implicit none
private
!Subroutines---------------------------------------------------------------
!Constructor
public :: ConstructRunOff
private :: AllocateInstance
private :: ReadDataFile
private :: AllocateVariables
private :: ConstructOverLandCoefficient
private :: ConstructStormWaterDrainage
private :: WriteStreetGutterLinksFile
private :: ConstructHDF5Output
private :: ConstructTimeSeries
!Selector
public :: GetOverLandFlow
public :: GetManning
public :: GetManningDelta
public :: GetFlowToChannels
public :: GetBoundaryImposed
public :: GetRouteDFour
public :: GetRouteDFourCells
public :: GetRouteDFourNeighbours
public :: GetRouteDFourFlux
public :: GetBoundaryFlux
public :: GetBoundaryCells
public :: GetFlowDischarge
public :: GetRunOffTotalDischargeFlowVolume
public :: GetRunoffWaterLevel
public :: GetRunoffWaterColumn !Final WaterColumn
public :: GetRunoffWaterColumnOld !Initial WaterColumn
public :: GetRunoffWaterColumnAT !WaterColumn After Transport (For RP)
public :: GetRunoffCenterVelocity
public :: GetRunoffTotalStoredVolume
public :: GetRunOffStoredVolumes
public :: GetRunOffBoundaryFlowVolume
public :: GetMassError
public :: GetNextRunOffDT
public :: SetBasinColumnToRunoff
public :: UnGetRunOff
!Modifier
public :: ModifyRunOff
private :: LocalWaterColumn
private :: IntegrateFlow
private :: ComputeNextDT
private :: RunOffOutput
private :: OutputTimeSeries
private :: AdjustSlope
!Destructor
public :: KillRunOff
!Management
private :: ReadLockExternalVar
private :: ReadUnLockExternalVar
private :: Ready
private :: LocateObjRunOff
!Interfaces----------------------------------------------------------------
private :: UnGetRunOff2D_R4
private :: UnGetRunOff2D_R8
interface UnGetRunOff
module procedure UnGetRunOff2D_R4
module procedure UnGetRunOff2D_R8
end interface UnGetRunOff
!Parameters----------------------------------------------------------------
integer, parameter :: KinematicWave_ = 1
integer, parameter :: DiffusionWave_ = 2
integer, parameter :: DynamicWave_ = 3
integer, parameter :: UnitMax = 80
!water column computation in faces
integer, parameter :: WCMaxBottom_ = 1
integer, parameter :: WCAverageBottom_ = 2
!Boundary flux
integer, parameter :: ComputeFlow_ = 1
integer, parameter :: InstantaneousFlow_ = 2
!Route D4 flux
integer, parameter :: Celerity_ = 1
integer, parameter :: Manning_ = 2
!Restart fiels format
integer, parameter :: BIN_ = 1
integer, parameter :: HDF_ = 2
!Types---------------------------------------------------------------------
!TODO: Use inheritance to make these types less code management and repititions
!GridPoint on top of river 1D node (will recieve water level from 1D model)
type T_NodeGridPoint
integer :: ID = null_int
integer :: GridI = null_int
integer :: GridJ = null_int
real :: RiverLevel = null_real
type(T_NodeGridPoint), pointer :: Next => null()
type(T_NodeGridPoint), pointer :: Prev => null()
end type
!GridPoint on top of river right and left banks (will recieve water level from associated NodeGridPoint
!and will be the cells where is computed river interaction flow)
type T_BankGridPoint
integer :: ID = null_int
integer :: GridI = null_int
integer :: GridJ = null_int
integer :: NGPId = null_int !associated NodeGridPoint
real :: RiverLevel = null_real
type(T_BankGridPoint), pointer :: Next => null()
type(T_BankGridPoint), pointer :: Prev => null()
end type
!GridPoint on margins (will recieve water level interpolated from associated BankGridPoint's
!and will integrate their river interaction flow into closest BGP - based on interpolation X)
type T_MarginGridPoint
integer :: ID = null_int
integer :: GridI = null_int
integer :: GridJ = null_int
integer :: BGPUpId = null_int !associated BankGridPoint upstream
integer :: BGPDownId = null_int !associated BankGridPoint downstream
real :: InterpolationFraction = null_real !x fraction from BGP upstream to downstream segment (at cell center)
integer :: BGPIntegrateFluxId = null_int !associated BGP where to route computed flow (BGP upstream or downstream)
integer :: NGPIntegrateFluxId = null_int !associated NGP where to route computed flow (NGP associated to BGP)
integer :: GridIIntegrateFlux = null_int !I where to integrate flux (BGP in case of DN or NGP in case OpenMI)
integer :: GridJIntegrateFlux = null_int !J where to integrate flux (BGP in case of DN or NGP in case OpenMI)
real :: RiverLevel = null_real
type(T_MarginGridPoint), pointer :: Next => null()
type(T_MarginGridPoint), pointer :: Prev => null()
end type
type T_OutPut
type (T_Time), pointer, dimension(:) :: OutTime => null()
integer :: NextOutPut = 1
logical :: Yes = .false.
type (T_Time), dimension(:), pointer :: RestartOutTime => null()
logical :: WriteRestartFile = .false.
logical :: RestartOverwrite = .false.
integer :: NextRestartOutput = 1
integer :: RestartFormat = BIN_
logical :: BoxFluxes = .false.
logical :: OutputFloodRisk = .false.
real :: FloodRiskVelCoef = null_real
logical :: WriteMaxFlowModulus = .false.
character(Pathlength) :: MaxFlowModulusFile = null_str
real, dimension(:,:), pointer :: MaxFlowModulus => null()
logical :: WriteMaxWaterColumn = .false.
character(Pathlength) :: MaxWaterColumnFile = null_str
real, dimension(:,:), pointer :: MaxWaterColumn => null()
logical :: WriteVelocityAtMaxWaterColumn = .false.
character(Pathlength) :: VelocityAtMaxWaterColumnFile = null_str
real, dimension(:,:), pointer :: VelocityAtMaxWaterColumn => null()
logical :: WriteMaxFloodRisk = .false.
character(Pathlength) :: MaxFloodRiskFile = null_str
real, dimension(:,:), pointer :: MaxFloodRisk => null()
logical :: WriteFloodPeriod = .false.
character(Pathlength) :: FloodPeriodFile = null_str
real, dimension(:,:), pointer :: FloodPeriod => null()
real :: FloodWaterColumnLimit = null_real
logical :: TimeSeries = .false.
logical :: TimeSerieDischON = .false.
integer :: DischargesNumber = null_int
integer, dimension(:), pointer :: TimeSerieDischID => null()
real, dimension(:,:), pointer :: TimeSerieDischProp => null()
integer :: TS_Numb_DischProp = null_int
type (T_Time) :: NextOutPutDisch
real :: OutPutDischDT
character(len=PathLength) :: TimeSerieLocationFile, DiscTimeSerieLocationFile
end type T_OutPut
type T_Files
character(PathLength) :: DataFile = null_str
character(PathLength) :: InitialFile = null_str
character(PathLength) :: FinalFile = null_str
character(PathLength) :: TransientHDF = null_str
character(PathLength) :: BoxesFile = null_str
end type T_Files
type T_ExtVar
integer, dimension(:,:), pointer :: BasinPoints => null()
real(8), dimension(:,:), pointer :: WaterColumn => null()
real , dimension(:,:), pointer :: GridCellArea => null()
real , dimension(:,:), pointer :: DUX, DUY => null()
real , dimension(:,:), pointer :: DVX, DVY => null()
real , dimension(:,:), pointer :: DXX, DYY => null()
real , dimension(:,:), pointer :: DZX, DZY => null()
real , dimension(:,:), pointer :: XX2D_Z, YY2D_Z => null()
real , dimension(:,:), pointer :: Topography => null()
integer, dimension(:,:), pointer :: BoundaryPoints2D => null()
integer, dimension(:,:), pointer :: RiverPoints => null()
real , dimension(:,:), pointer :: CellSlope => null()
type (T_Time) :: Now
real :: DT = null_real
end type T_ExtVar
type T_Converge
integer :: MinIterations = 1
integer :: MaxIterations = 1024
logical :: Stabilize = .false.
real :: StabilizeFactor = 0.01
real :: DTFactorUp = 1.25
real :: DTFactorDown = 1.25
real :: StabilizeHardCutLimit = 128
real :: DTSplitFactor = 2.0
real :: CurrentDT = null_real
real :: NextDT = null_real
integer :: LastGoodNiteration = 1
integer :: NextNiteration = 1
logical :: LimitDTCourant = .false.
real :: MaxCourant = 1.0
integer :: MinToRestart = 0
real :: MinimumValueToStabilize = 0.001
logical :: CheckDecreaseOnly = .false.
end type T_Converge
type T_FromTimeSerie
integer :: ObjTimeSerie = 0
character(len=PathLength) :: FileName = null_str
integer :: DataColumn = null_int
end type T_FromTimeSerie
!level imposed as time serie
type T_ImposedLevelTS
character(len=StringLength) :: Name = null_str
! type(T_PointF) :: Location
integer :: ID = null_int
real :: DefaultValue = null_real
character(len=StringLength) :: ValueType = null_str
logical :: TimeSerieHasData = .false.
type(T_FromTimeSerie) :: TimeSerie
end type T_ImposedLevelTS
type T_RunOff
integer :: InstanceID = 0
character(len=StringLength) :: ModelName = null_str
integer :: ObjBasinGeometry = 0
integer :: ObjTime = 0
integer :: ObjHorizontalGrid = 0
integer :: ObjHorizontalMap = 0
integer :: ObjGridData = 0
integer :: ObjHDF5 = 0
integer :: ObjDrainageNetwork = 0
integer :: ObjDischarges = 0
integer :: ObjEnterData = 0
integer :: ObjBoxDif = 0
integer :: ObjTimeSerie = 0
type (T_OutPut ) :: OutPut
type (T_ExtVar) :: ExtVar
type (T_Files) :: Files
type (T_Time) :: BeginTime
type (T_Time) :: EndTime
real(8), dimension(:,:), pointer :: myWaterLevel => null()
real(8), dimension(:,:), pointer :: myWaterColumn => null()
real, dimension(:,:), pointer :: InitialWaterColumn => null()
real, dimension(:,:), pointer :: InitialWaterLevel => null()
logical :: PresentInitialWaterColumn = .false.
logical :: PresentInitialWaterLevel = .false.
real(8), dimension(:,:), pointer :: myWaterVolume => null()
real(8), dimension(:,:), pointer :: myWaterColumnOld => null() !OldColumn from Basin
real(8), dimension(:,:), pointer :: myWaterColumnAfterTransport => null() !for property transport
real(8), dimension(:,:), pointer :: myWaterVolumePred => null() !to avoid negative collumns
real(8), dimension(:,:), pointer :: myWaterVolumeOld => null()
real, dimension(:,:), pointer :: lFlowToChannels => null() !Instantaneous Flow To Channels
real, dimension(:,:), pointer :: iFlowToChannels => null() !Integrated Flow
real, dimension(:,:), pointer :: iFlowRouteDFour => null() !Integrated Route D4 flux
real, dimension(:,:), pointer :: lFlowBoundary => null() !Instantaneous Flow to impose BC
real, dimension(:,:), pointer :: iFlowBoundary => null() !Integrated Flow to impose BC
real, dimension(:,:), pointer :: lFlowDischarge => null() !Instantaneous Flow of discharges
real, dimension(:,:), pointer :: iFlowDischarge => null() !Integrated Flow of discharges
real(8), dimension(:,:), pointer :: lFlowX, lFlowY => null() !Instantaneous OverLandFlow (LocalDT )
real(8), dimension(:,:), pointer :: iFlowX, iFlowY => null() !Integrated OverLandFlow (AfterSumDT)
real(8), dimension(:,:), pointer :: FlowXOld, FlowYOld => null() !Flow From previous iteration
real(8), dimension(:,:), pointer :: InitialFlowX, InitialFlowY => null() !Initial Flow of convergence
real(8), dimension(:,:), pointer :: VelModFaceU, VelModFaceV => null() !Flow From previous iteration
real, dimension(:,:), pointer :: AreaU, AreaV => null()
integer, dimension(:,:), pointer :: ComputeFaceU => null()
integer, dimension(:,:), pointer :: ComputeFaceV => null()
integer, dimension(:,:), pointer :: ComputeAdvectionU => null()
integer, dimension(:,:), pointer :: ComputeAdvectionV => null()
real, dimension(:,:), pointer :: NoAdvectionPoints => null()
integer, dimension(:,:), pointer :: OpenPoints => null() !Mask for gridcells above min watercolumn
real, dimension(:,:), pointer :: OverLandCoefficient => null() !Manning or Chezy
real, dimension(:,:), pointer :: OverLandCoefficientDelta => null() !For erosion/deposition
real, dimension(:,:), pointer :: OverLandCoefficientX => null() !Manning or Chezy
real, dimension(:,:), pointer :: OverLandCoefficientY => null() !Manning or Chezy
real, dimension(:,:), pointer :: StormWaterDrainageCoef => null() !Sewer System Percentagem (area)
real, dimension(:,:), pointer :: StormWaterVolume => null() !Volume of storm water stored in each
!cell
real, dimension(:,:), pointer :: StormWaterFlowX => null() !Auxilizary Var for explicit routing
real, dimension(:,:), pointer :: StormWaterFlowY => null() !Auxilizary Var for explicit routing
real, dimension(:,:), pointer :: StormWaterCenterFlowX => null() !Output
real, dimension(:,:), pointer :: StormWaterCenterFlowY => null() !Output
real, dimension(:,:), pointer :: StormWaterCenterModulus => null() !Output
real, dimension(:,:), pointer :: BuildingsHeight => null() !Height of building in cell
real, dimension(:,:), pointer :: NumberOfSewerStormWaterNodes => null() !Number of total SWMM nodes
!(sewer + storm water) per grid cell
!that interact with MOHID
real, dimension(:,:), pointer :: NumberOfStormWaterNodes => null() !Number of SWMM storm water only nodes
!per grid cell that interact
!with MOHID (default is the same as
!NumberOfSewerStormWaterNodes)
real, dimension(:,:), pointer :: StreetGutterLength => null() !Length of Stret Gutter in a given cell
real, dimension(:,:), pointer :: MassError => null() !Contains mass error
real, dimension(:,:), pointer :: CenterFlowX => null()
real, dimension(:,:), pointer :: CenterFlowY => null()
real, dimension(:,:), pointer :: CenterVelocityX => null()
real, dimension(:,:), pointer :: CenterVelocityY => null()
real, dimension(:,:), pointer :: FlowModulus => null()
real, dimension(:,:), pointer :: VelocityModulus => null()
integer, dimension(:,:), pointer :: LowestNeighborI => null() !Lowest Neighbor in the surroundings
integer, dimension(:,:), pointer :: LowestNeighborJ => null() !Lowest Neighbor in the surroundings
integer, dimension(:,:), pointer :: DFourSinkPoint => null() !Point which can't drain with in X/Y only
integer, dimension(:,:), pointer :: StabilityPoints => null() !Points where models check stability
type(T_PropertyID) :: OverLandCoefficientID, NoAdvectionZonesID
logical :: StormWaterModel = .false. !If connected to SWMM
real, dimension(:,:), pointer :: StormWaterEffectiveFlow => null() !Flow from SWMM (inflow + outflow)
real, dimension(:,:), pointer :: StreetGutterPotentialFlow=> null() !Potential flow to street gutters
!in grid cells with street gutters
real, dimension(:,:), pointer :: StormWaterPotentialFlow => null() !Sum of all potential flows
!from street gutters draining to
!grid cells with storm water nodes
real, dimension(:,:), pointer :: StreetGutterEffectiveFlow=> null() !Effective flow to street gutters
!in grid cells with street gutters
integer, dimension(:,:), pointer :: StreetGutterTargetI => null() !Sewer interaction point...
integer, dimension(:,:), pointer :: StreetGutterTargetJ => null() !...where street gutter drains to
real, dimension(:,:), pointer :: NodeRiverLevel => null() !river level at river points (from DN or external model)
integer, dimension(:,:), pointer :: NodeRiverMapping => null() !mapping of river points where interaction occurs (for external model)
real :: MinSlope = null_real
logical :: AdjustSlope = .false.
logical :: Stabilize = .false.
logical :: Discharges = .false.
logical :: RouteDFourPoints = .false.
logical :: RouteDFourPointsOnDN = .false.
integer :: RouteDFourMethod = null_int
logical :: StormWaterDrainage = .false.
real :: StormWaterInfiltrationVelocity = 1.4e-5 !~50mm/h
real :: StormWaterFlowVelocity = 0.2 !velocity in pipes
logical :: Buildings = .false.
! real :: StabilizeFactor = null_real
! real :: StabilizeHardCutLimit = 0.1
integer :: HydrodynamicApproximation = DiffusionWave_
logical :: CalculateAdvection = .true.
logical :: NoAdvectionZones = .false.
logical :: CalculateDiffusion = .false.
real :: Viscosity_Coef = null_real
logical :: CalculateCellMargins = .true.
logical :: ImposeMaxVelocity = .false.
real :: ImposedMaxVelocity = 0.1
! integer :: LastGoodNiter = 1
! integer :: NextNiter = 1
! real :: InternalTimeStepSplit = 1.5
! integer :: MinIterations = 1
! integer :: MinToRestart = 0
real :: MinimumWaterColumn = null_real
real :: MinimumWaterColumnAdvection = null_real
! real :: MinimumWaterColumnStabilize = null_real
! real :: NextDT = null_real
! real :: DTFactor = null_real
! real :: DTFactorUp = null_real
! real :: DTFactorDown = null_real
! real :: CurrentDT = null_real
! logical :: LimitDTCourant = .false.
! logical :: LimitDTVariation = .true.
! real :: MaxCourant = 1.0
logical :: ImposeBoundaryValue = .false.
logical :: AllowBoundaryInflow = .false.
logical :: BoundaryImposedLevelInTime = .false.
real :: BoundaryValue = null_real
real :: MaxDtmForBoundary = null_real
integer :: BoundaryMethod = null_int
integer, dimension(:,:), pointer :: BoundaryCells => null()
type (T_ImposedLevelTS) :: ImposedLevelTS
! integer :: MaxIterations = 5
logical :: SimpleChannelInteraction = .false.
logical :: LimitToCriticalFlow = .true.
integer :: FaceWaterColumn = WCMaxBottom_
! real :: MaxVariation = null_real
integer :: OverlandChannelInteractionMethod = null_int
! logical :: CheckDecreaseOnly = .false.
type(T_Converge) :: CV !Convergence options
real(8) :: BoundaryFlowVolume = 0.0 !m3 => positive if flow is towards boundary.
real(8) :: VolumeStoredInSurface = 0.0
real(8) :: VolumeStoredInStormSystem = 0.0
real(8) :: TotalDischargeFlowVolume = 0.0
logical :: Continuous = .false.
logical :: StopOnWrongDate = .true.
real(8) :: TotalStoredVolume = 0.
integer :: BasinCellsCount = 0
!Grid size
type (T_Size2D) :: Size
type (T_Size2D) :: WorkSize
type(T_NodeGridPoint ), pointer :: FirstNodeGridPoint => null()
type(T_NodeGridPoint ), pointer :: LastNodeGridPoint => null()
integer :: NodeGridPointNumber = 0
type(T_MarginGridPoint ), pointer :: FirstMarginGridPoint => null()
type(T_MarginGridPoint ), pointer :: LastMarginGridPoint => null()
integer :: MarginGridPointNumber = 0
type(T_BankGridPoint ), pointer :: FirstBankGridPoint => null()
type(T_BankGridPoint ), pointer :: LastBankGridPoint => null()
integer :: BankGridPointNumber = 0
logical :: Use1D2DInteractionMapping = .false.
type(T_RunOff), pointer :: Next => null()
end type T_RunOff
!Global Module Variables
type (T_RunOff), pointer :: FirstObjRunOff => null()
type (T_RunOff), pointer :: Me => null()
!--------------------------------------------------------------------------
contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine ConstructRunOff(ModelName, &
RunOffID, &
ComputeTimeID, &
HorizontalGridID, &
HorizontalMapID, &
GridDataID, &
BasinGeometryID, &
DrainageNetworkID, &
DischargesID, &
STAT)
!Arguments---------------------------------------------------------------
character(len=*) :: ModelName
integer :: RunOffID
integer :: ComputeTimeID
integer :: HorizontalGridID
integer :: HorizontalMapID
integer :: GridDataID
integer :: BasinGeometryID
integer :: DrainageNetworkID
integer, optional, intent(OUT) :: STAT
integer, intent (OUT) :: DischargesID
!External----------------------------------------------------------------
integer :: ready_
!Local-------------------------------------------------------------------
integer :: STAT_, STAT_CALL
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mRunOff_)) then
nullify (FirstObjRunOff)
call RegisterModule (mRunOff_)
endif
call Ready(RunOffID, ready_)
cd0 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
Me%ModelName = ModelName
!Associates External Instances
Me%ObjTime = AssociateInstance (mTIME_ , ComputeTimeID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_ , HorizontalGridID )
Me%ObjHorizontalMap = AssociateInstance (mHORIZONTALMAP_ , HorizontalMapID )
Me%ObjGridData = AssociateInstance (mGRIDDATA_ , GridDataID )
Me%ObjBasinGeometry = AssociateInstance (mBASINGEOMETRY_ , BasinGeometryID )
if (DrainageNetworkID /= 0) then
Me%ObjDrainageNetwork = AssociateInstance (mDRAINAGENETWORK_, DrainageNetworkID)
endif
!Time Stuff
call GetComputeTimeLimits (Me%ObjTime, BeginTime = Me%BeginTime, &
EndTime = Me%EndTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructRunOff - ModuleRunOff - ERR010'
call GetComputeTimeStep (Me%ObjTime, Me%ExtVar%DT, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructRunOff - ModuleRunOff - ERR011'
Me%CV%NextNiteration = 1
Me%CV%CurrentDT = Me%ExtVar%DT
call ReadLockExternalVar (StaticOnly = .false.)
!Gets the size of the grid
call GetHorizontalGridSize (Me%ObjHorizontalGrid, &
Size = Me%Size, &
WorkSize = Me%WorkSize, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructRunOff - ModuleRunOff - ERR020'
call AllocateVariables
call ReadDataFile
call InitializeVariables
call ConstructOverLandCoefficient
!Checks if River Network is consistent with the one previously constructed
if (DrainageNetworkID /= 0) then
call CheckRiverNetWorkConsistency
endif
!Constructs Discharges
if (Me%Discharges) then
call ConstructDischarges
endif
!Constructs StormWaterDrainage
if (Me%StormWaterDrainage .or. Me%StormWaterModel) then
call ConstructStormWaterDrainage
endif
!Constructs Boundary Cells
if (Me%ImposeBoundaryValue) then
call CheckBoundaryCells
if (Me%BoundaryImposedLevelInTime) call ModifyBoundaryLevel
endif
!Reads conditions from previous run
if (Me%Continuous) then
if (Me%OutPut%RestartFormat == BIN_) then
call ReadInitialFile_Bin
else if (Me%OutPut%RestartFormat == HDF_) then
call ReadInitialFile_Hdf
endif
endif
if (Me%OutPut%Yes) then
call ConstructHDF5Output
endif
call CalculateTotalStoredVolume
!Output Results
if (Me%OutPut%Yes .or. Me%OutPut%TimeSeries) then
call ComputeCenterValues
endif
if(Me%OutPut%Yes)then
call RunOffOutput
endif
if(Me%OutPut%TimeSeries) then
call OutputTimeSeries
endif
call ReadUnLockExternalVar (StaticOnly = .false.)
!Returns ID
RunOffID = Me%InstanceID
DischargesID = Me%ObjDischarges
STAT_ = SUCCESS_
else cd0
stop 'ModuleRunOff - ConstructRunOff - ERR030'
end if cd0
if (present(STAT)) STAT = STAT_
!----------------------------------------------------------------------
end subroutine ConstructRunOff
!--------------------------------------------------------------------------
subroutine AllocateInstance
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
type (T_RunOff), pointer :: NewObjRunOff
type (T_RunOff), pointer :: PreviousObjRunOff
!Allocates new instance
allocate (NewObjRunOff)
nullify (NewObjRunOff%Next)
!Insert New Instance into list and makes Current point to it
if (.not. associated(FirstObjRunOff)) then
FirstObjRunOff => NewObjRunOff
Me => NewObjRunOff
else
PreviousObjRunOff => FirstObjRunOff
Me => FirstObjRunOff%Next
do while (associated(Me))
PreviousObjRunOff => Me
Me => Me%Next
enddo
Me => NewObjRunOff
PreviousObjRunOff%Next => NewObjRunOff
endif
Me%InstanceID = RegisterNewInstance (mRUNOFF_)
end subroutine AllocateInstance
!--------------------------------------------------------------------------
subroutine ReadDataFile
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: STAT_CALL
type(T_PropertyID) :: InitialWaterColumnID, InitialWaterLevelID
type(T_PropertyID) :: OverLandCoefficientDeltaID
type(T_PropertyID) :: StormWaterDrainageID
type(T_PropertyID) :: BuildingsHeightID
type(T_PropertyID) :: NumberOfSewerStormWaterNodesID
type(T_PropertyID) :: NumberOfStormWaterNodesID
type(T_PropertyID) :: StreetGutterLengthID
integer :: ObjEnterDataGutterInteraction = 0
character(len=StringLength) :: InitializationMethod
character(len=PathLength) :: Filename, MappingFileName
character(len=StringLength) :: StormWaterGutterRegExpression, StormWaterGutterRegExpressionFromGD
integer :: iflag, ClientNumber, FoundSWMMRegExpression
logical :: BlockFound
integer :: i, j
logical :: DynamicAdjustManning
real :: dummy
!Reads the name of the data file from nomfich
call ReadFileName ('RUNOFF_DATA', Me%Files%DataFile, "RunOff Data File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR010'
!Reads the name of the transient HDF file from nomfich
call ReadFileName ('RUNOFF_HDF', Me%Files%TransientHDF, "RunOff HDF File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR020'
call ReadFileName('RUNOFF_FIN', Me%Files%FinalFile, &
Message = "RunOff Final File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR030'
!Constructs the DataFile
call ConstructEnterData (Me%ObjEnterData, Me%Files%DataFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR040'
!Initial Water Column
call GetData(dummy, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'INITIAL_WATER_COLUMN', &
default = 0.0, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR050'
if (iflag /= 0) then
write(*,*)'The keyword INITIAL_WATER_COLUMN is obsolete.'
write(*,*)'Please use the block <BeginInitialWaterColumn> / <EndInitialWaterColumn>'
stop 'ReadDataFile - ModuleRunOff - ERR060'
endif
!Gets Block
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginInitialWaterColumn>', &
'<EndInitialWaterColumn>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR070'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = InitialWaterColumnID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%InitialWaterColumn, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR080'
call KillFillMatrix(InitialWaterColumnID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR081'
Me%PresentInitialWaterColumn = .true.
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR090'
else
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginInitialWaterLevel>', &
'<EndInitialWaterLevel>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR091'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = InitialWaterLevelID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%InitialWaterLevel, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR092'
call KillFillMatrix(InitialWaterLevelID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR093'
Me%PresentInitialWaterLevel = .true.
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR094'
else
write(*,*)
write(*,*)'Missing Block <BeginInitialWaterColumn> / <EndInitialWaterColumn>'
write(*,*)'or <BeginInitialWaterLevel> / <EndInitialWaterLevel>'
stop 'ReadDataFile - ModuleRunOff - ERR100'
endif
endif
!Gets Minimum Slope
call GetData(Me%MinSlope, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'MIN_SLOPE', &
default = 0.0, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0110'
if (Me%MinSlope < 0.0 .or. Me%MinSlope >= 1.) then
write (*,*) 'Invalid Minimum Slope [MIN_SLOPE]'
stop 'ReadDataFile - ModuleRunOff - ERR0120'
end if
!Adjusts Slope according to
!http://www.hkh-friend.net.np/rhdc/training/lectures/HEGGEN/Tc_3.pdf
call GetData(Me%AdjustSlope, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'ADJUST_SLOPE', &
default = .true., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0130'
!Gets Routing method
call GetData(Me%HydrodynamicApproximation, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'HYDRODYNAMIC_APROX', &
default = DiffusionWave_, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0140'
if (Me%HydrodynamicApproximation /= KinematicWave_ .and. &
Me%HydrodynamicApproximation /= DiffusionWave_ .and. &
Me%HydrodynamicApproximation /= DynamicWave_) then
write (*,*) 'Invalid Hydrodynamic Approximation [HYDRODYNAMIC_APROX]'
stop 'ReadDataFile - ModuleRunOff - ERR0150'
end if
if (Me%HydrodynamicApproximation == DynamicWave_) then
!Gets if advection is to be calculated
call GetData(Me%CalculateAdvection, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'ADVECTION', &
default = .true., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0160'
if (Me%CalculateAdvection) then
!Minimum Water Column for advection computation
call GetData(Me%MinimumWaterColumnAdvection, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'MIN_WATER_COLUMN_ADVECTION', &
default = 0.0, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0170'
call GetData(Me%NoAdvectionZones, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'NO_ADVECTION_ZONES', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0170a'
if(Me%NoAdvectionZones)then
call ConstructAdvectionZones
else
call SetMatrixValue(Me%ComputeAdvectionU, Me%Size, 1)
call SetMatrixValue(Me%ComputeAdvectionU, Me%Size, 1)
endif
else
call SetMatrixValue(Me%ComputeAdvectionU, Me%Size, 0)
call SetMatrixValue(Me%ComputeAdvectionU, Me%Size, 0)
endif
endif
!Method for computing water column in the face (1 - Using max level and max bottom;
!2- using max level and average of bottom)
call GetData(Me%FaceWaterColumn, &
Me%ObjEnterData, iflag, &
keyword = 'WATER_COLUMN_FACE', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = WCMaxBottom_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR175'
if (Me%FaceWaterColumn /= WCMaxBottom_ .and. Me%FaceWaterColumn /= WCAverageBottom_) then
write(*,*) 'Unknown option for WATER_COLUMN_FACE'
stop 'ReadDataFile - ModuleRunOff - ERR176'
endif
if (Me%FaceWaterColumn == WCMaxBottom_) then
!Gets if compute "margins" aside of adjacent cells that produce friction
call GetData(Me%CalculateCellMargins, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'HYDRAULIC_RADIUS_MARGINS', &
default = .true., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR180'
endif
!Gets if solution is limited by an maximum velocity
call GetData(Me%ImposeMaxVelocity, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'IMPOSE_MAX_VELOCITY', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR190'
if (Me%ImposeMaxVelocity) then
!Gets if solution is limited by an maximum velocity
call GetData(Me%ImposedMaxVelocity, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'MAX_VELOCITY', &
default = 0.1, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR200'
endif
!Gets if Manning Coeficient is increased with water depth
call GetData(DynamicAdjustManning, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'DYNAMIC_ADJUST_MANNING', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR210'
if (iflag > 0 .and. .not. DynamicAdjustManning) then
write(*,*)'The option DynamicAdjustManning (DYNAMIC_ADJUST_MANNING) has been removed.'
write(*,*)'Please review your runoff data file!'
endif
if (DynamicAdjustManning) then
write(*,*)'The option DynamicAdjustManning (DYNAMIC_ADJUST_MANNING) has been removed.'
write(*,*)'Please review your runoff data file!'
stop
endif
!Minimum Water Column for overland flow
call GetData(Me%MinimumWaterColumn, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'MIN_WATER_COLUMN', &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR220'
if (iflag == 0) then
write(*,*)'MIN_WATER_COLUMN must be defined in module Runoff'
stop 'ReadDataFile - ModuleRunOff - ERR230'
endif
!Continuous Computation
call GetData(Me%Continuous, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'CONTINUOUS', &
default = .false., &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR240'
if (Me%Continuous) then
call ReadFileName('RUNOFF_INI', Me%Files%InitialFile, &
Message = "Runoff Initial File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0250'
endif
call GetData(Me%StopOnWrongDate, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'STOP_ON_WRONG_DATE', &
default = .true., &
ClientModule = 'ModuleBasin', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR260'
!Impose Boundary Value
call GetData(Me%ImposeBoundaryValue, &
Me%ObjEnterData, iflag, &
keyword = 'IMPOSE_BOUNDARY_VALUE', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR350'
if (Me%ImposeBoundaryValue) then
!verify if the user wants to allow water to go in the domain (true) if
!boundary level higher than water level or not (false) and the level imposed
!behaves like a wall, only exits if higher and does not allow to get inside
call GetData(Me%AllowBoundaryInflow, &
Me%ObjEnterData, iflag, &
keyword = 'ALLOW_BOUNDARY_INFLOW', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR363'
call ReadBoundaryConditions
endif
!Discharges
call GetData(Me%Discharges, &
Me%ObjEnterData, iflag, &
keyword = 'DISCHARGES', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR370'
!Discharges output time series
if (Me%Discharges) then
call GetData(Me%Output%TimeSerieDischON, &
Me%ObjEnterData, iflag, &
keyword = 'TIME_SERIE_DISCHARGES', &
Default = .false., &
SearchType = FromFile, &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR375'
if (Me%Output%TimeSerieDischON) then
call GetData(Me%Output%DiscTimeSerieLocationFile, &
Me%ObjEnterData,iflag, &
SearchType = FromFile, &
keyword = 'DISCHARGE_TIME_SERIE_LOCATION', &
ClientModule = 'ModuleRunOff', &
Default = Me%Files%DataFile, &
STAT = STAT_CALL)
if (STAT_CALL/=SUCCESS_)stop 'ReadDataFile - ModuleRunOff - ERR377'
Me%Output%NextOutPutDisch = Me%BeginTime
call GetData(Me%Output%OutPutDischDT, &
Me%ObjEnterData, iflag, &
keyword = 'DISCHARGE_DT_OUTPUT_TIME', &
Default = FillValueReal, &
SearchType = FromFile, &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL/=SUCCESS_)stop 'ReadDataFile - ModuleRunOff - ERR378'
if (iflag == 0) then
call GetComputeTimeStep(Me%ObjTime, &
Me%Output%OutPutDischDT, &
STAT = STAT_CALL)
if (STAT_CALL/=SUCCESS_)stop 'ReadDataFile - ModuleRunOff - ERR379'
endif
endif
endif
!Discharges
call GetData(Me%SimpleChannelInteraction, &
Me%ObjEnterData, iflag, &
keyword = 'SIMPLE_CHANNEL_FLOW', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR380'
!Routes D4 Points
call GetData(Me%RouteDFourPoints, &
Me%ObjEnterData, iflag, &
keyword = 'ROUTE_D4', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR390'
if (Me%RouteDFourPoints) then
!Routes D4 Points
call GetData(Me%RouteDFourPointsOnDN, &
Me%ObjEnterData, iflag, &
keyword = 'ROUTE_D4_ON_DN', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR391'
call GetData(Me%RouteDFourMethod, &
Me%ObjEnterData, iflag, &
keyword = 'ROUTE_D4_METHOD', &
! Default = Celerity_, &
Default = Manning_, &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR392'
if (Me%RouteDFourMethod /= Celerity_ .and. Me%RouteDFourMethod /= Manning_) then
write(*,*)'ROUTE_D4_METHOD must be or 1 - Celerity based or 2 - Manning Equation'
stop 'ReadDataFile - ModuleRunOff - ERR0393'
endif
endif
!Limits Flow to critical
call GetData(Me%LimitToCriticalFlow, &
Me%ObjEnterData, iflag, &
keyword = 'LIMIT_TO_CRITICAL_FLOW', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .true., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR394'
!Storm Water Drainage
call GetData(Me%StormWaterDrainage, &
Me%ObjEnterData, iflag, &
keyword = 'STORM_WATER', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR400'
if (Me%StormWaterDrainage) then
!Storm Water Infiltration Velocity
call GetData(Me%StormWaterInfiltrationVelocity, &
Me%ObjEnterData, iflag, &
keyword = 'STORM_WATER_INF_VEL', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
default = 1.4e-5, & !~50mm/h
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR410'
!Storm Water Transfer Coeficient
call GetData(Me%StormWaterFlowVelocity, &
Me%ObjEnterData, iflag, &
keyword = 'STORM_WATER_FLOW_VEL', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
default = 0.2, & !0.2m/s
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR420'
endif
!If Buildings are to be simulated (flow ocuation in urban areas)
call GetData(Me%Buildings, &
Me%ObjEnterData, iflag, &
keyword = 'BUILDINGS', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR430'
!If Connected to a StormWater model
call GetData(Me%StormWaterModel, &
Me%ObjEnterData, iflag, &
keyword = 'STORM_WATER_MODEL_LINK', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR440'
!Gets Output Time
call GetOutPutTime(Me%ObjEnterData, &
CurrentTime = Me%ExtVar%Now, &
EndTime = Me%EndTime, &
keyword = 'OUTPUT_TIME', &
SearchType = FromFile, &
OutPutsTime = Me%OutPut%OutTime, &
OutPutsOn = Me%OutPut%Yes, &
STAT = STAT_CALL)
Me%OutPut%NextOutPut = 1
!Output for restart
call GetOutPutTime(Me%ObjEnterData, &
CurrentTime = Me%ExtVar%Now, &
EndTime = Me%EndTime, &
keyword = 'RESTART_FILE_OUTPUT_TIME', &
SearchType = FromFile, &
OutPutsTime = Me%OutPut%RestartOutTime, &
OutPutsOn = Me%OutPut%WriteRestartFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunoff - ERR450'
call GetData(Me%OutPut%RestartFormat, &
Me%ObjEnterData, &
iflag, &
SearchType = FromFile, &
keyword = 'RESTART_FILE_FORMAT', &
Default = BIN_, &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunoff - ERR452'
if (Me%OutPut%RestartFormat /= BIN_ .and. Me%OutPut%RestartFormat /= HDF_) then
write (*,*)
write (*,*) 'RESTART_FILE_FORMAT options are: 1 - Binary or 2 - HDF'
stop 'ReadDataFile - ModuleRunoff - ERR455'
endif
call GetData(Me%OutPut%RestartOverwrite, &
Me%ObjEnterData, &
iflag, &
SearchType = FromFile, &
keyword = 'RESTART_FILE_OVERWRITE', &
Default = .true., &
ClientModule = 'ModuleBasin', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunoff - ERR460'
call RewindBuffer (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0470'
!Gets Block for OverLand Coef
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginOverLandCoefficient>', &
'<EndOverLandCoefficient>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0480'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = Me%OverLandCoefficientID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%OverLandCoefficient, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0490'
call KillFillMatrix(Me%OverLandCoefficientID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0500'
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0501'
!Check that manning values entered are not zero or negative
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
if (.not. Me%OverLandCoefficient(i,j) .gt. 0.0) then
write(*,*) 'Found Manning Overland coefficient zero or negative in input'
write(*,*) 'in cell', i, j
stop 'ReadDataFile - ModuleRunoff - ERR0510'
endif
endif
enddo
enddo
else
write(*,*)'Missing Block <BeginOverLandCoefficient> / <EndOverLandCoefficient>'
stop 'ReadDataFile - ModuleRunOff - ERR0520'
endif
!Gets Block for OverLand Coef Difference
!To compute overland resistance in bottom for shear computation (erosion/deposition).
!This process was created to remove from manning the resistance given by
!aerial vegetation parts that affect flow but do not affect bottom shear. Without that,
!a manning increase (e.g. forestation) in one cell increases water depth (and reduces velocity)
!but may increase shear stress (because water height increase is transformed in bottom resistance
!using manning - chezy see module runoff properties)
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginOverLandCoefficientDelta>', &
'<EndOverLandCoefficientDelta>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0530'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = OverLandCoefficientDeltaID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%OverLandCoefficientDelta, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0540'
call KillFillMatrix(OverLandCoefficientDeltaID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0550'
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0551'
!Check that final manning values are not zero or negative
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
if (.not. (Me%OverLandCoefficient(i,j) - Me%OverLandCoefficientDelta(i,j)) .gt. 0.0) then
write(*,*) 'Manning Overland coefficient delta found zero or negative in input'
write(*,*) 'in cell', i, j
stop 'ReadDataFile - ModuleRunoff - ERR0560'
endif
endif
enddo
enddo
else
!Do not remove aerial vegetation effect from manning
Me%OverLandCoefficientDelta(:,:) = 0.0
endif
!Looks for StormWater DrainageCoef
if (Me%StormWaterDrainage) then
allocate(Me%StormWaterDrainageCoef (Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB))
Me%StormWaterDrainageCoef = null_real
call RewindBuffer (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR570'
!Gets Flag with Sewer Points
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginStormWaterDrainage>', &
'<EndStormWaterDrainage>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR580'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = StormWaterDrainageID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%StormWaterDrainageCoef, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR590'
call KillFillMatrix(StormWaterDrainageID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR600'
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0601'
else
write(*,*)'Missing Block <BeginStormWaterDrainage> / <EndStormWaterDrainage>'
stop 'ReadDataFile - ModuleRunOff - ERR0610'
endif
endif
allocate(Me%BuildingsHeight(Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB))
Me%BuildingsHeight = 0.0
if (Me%Buildings) then
call RewindBuffer (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR620'
!Gets Flag with Sewer Points
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginBuildingsHeight>', &
'<EndBuildingsHeight>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR630'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = BuildingsHeightID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%BuildingsHeight, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR640'
call KillFillMatrix(BuildingsHeightID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR650'
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0651'
else
write(*,*)'Missing Block <BeginBuildingsHeight> / <EndBuildingsHeight>'
stop 'ReadDataFile - ModuleRunOff - ERR0670'
endif
endif
!Looks for StormWater Interaction Point
if (Me%StormWaterModel) then
allocate(Me%NumberOfSewerStormWaterNodes(Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB))
allocate(Me%StreetGutterLength (Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB))
allocate(Me%NumberOfStormWaterNodes (Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB))
call RewindBuffer (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR680'
!Gets Flag with Sewer Points
!These are all cells that interact with SWMM (via SWMM junction (manhole) overflow to MOHID Land
!or MOHID Land street gutter inflow to SWMM - directed to nearest manhole)
!The cell value is number of manholes in each cell
!If this field is zero everywhere, there wil be no SWMM-MOHIDLand interaction
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginStormWaterInteraction>', &
'<EndStormWaterInteraction>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR690'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = NumberOfSewerStormWaterNodesID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%NumberOfSewerStormWaterNodes, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR691'
call KillFillMatrix(NumberOfSewerStormWaterNodesID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR692'
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0693'
else
write(*,*)'Missing Block <BeginStormWaterInteraction> / <EndStormWaterInteraction>'
stop 'ReadDataFile - ModuleRunOff - ERR694'
endif
call RewindBuffer (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR695'
!Look for keyword to filter SWMM nodes (regular expresion on node names) that can recieve MOHID Land flow (e.g. pluvial
!nodes) This can only be used in conjunction with <BeginStormWaterGutterInteraction>/<EndStormWaterGutterInteraction>
!to have gutter interaction number of nodes. this latter griddata needs to have been built with same regular exression
!(saved in grid data comment2 during grid data build with tool)
call GetData(StormWaterGutterRegExpression, &
Me%ObjEnterData, &
FoundSWMMRegExpression, &
SearchType = FromFile, &
keyword = 'SWMM_JUNCTIONS_GUTTER_REGEX', &
ClientModule = 'ModuleBasin', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunoff - ERR694.5'
!Gets Sewer Points that can interact with street gutter
!By default these are the same points as NumberOfSewerStormWaterNodes (all points can recieve street gutter)
!This exists to individualize SWMM junctions (manholes) that can recieve gutter flow (eg. pluvial junctions)
!It is only used to find for each gutter the closer cell that can have gutter inflow (avoid the ones that can't)
!If this field is zero everywhere, it cant find SWMM junctions to discharge gutter flow and will return error
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginStormWaterGutterInteraction>', &
'<EndStormWaterGutterInteraction>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR694'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = NumberOfStormWaterNodesID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%NumberOfStormWaterNodes, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFi le - ModuleRunOff - ERR695'
!verify that griddata file was created with same regular expression that SWMM will apply to
if (FoundSWMMRegExpression == 0) then
write(*,*)'With <BeginStormWaterGutterInteraction>/<EndStormWaterGutterInteraction>'
write(*,*)'SWMM_JUNCTIONS_GUTTER_REGEX must be defined in module Runoff.'
stop 'ReadDataFile - ModuleRunOff - ERR0694.6'
else
!open grid data and check what regular expression was used
call GetData(InitializationMethod, &
Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'INITIALIZATION_METHOD', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'ReadDataFile - ModuleRunOff - ERR694.7'
select case (trim(adjustl(InitializationMethod)))
case ("ASCII_File", "ASCII_FILE", "ascii_file", "Ascii_file")
call GetData(FileName, &
Me%ObjEnterData,iflag, &
SearchType = FromBlock, &
keyword = 'FILENAME', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'ReadDataFile - ModuleRunOff - ERR694.8'
call ConstructEnterData (ObjEnterDataGutterInteraction, FileName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR694.9'
!in comment 2 is grid data regular expression used to build the grid data
!(method that builds the grid data automatically writes that)
call GetData(StormWaterGutterRegExpressionFromGD, &
ObjEnterDataGutterInteraction,iflag, &
SearchType = FromFile, &
keyword = 'COMENT2', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'ReadDataFile - ModuleRunOff - ERR694.9'
if (StormWaterGutterRegExpression /= StormWaterGutterRegExpressionFromGD) then
write(*,*) 'Grid Data in <BeginStormWaterGutterInteraction>/<EndStormWaterGutterInteraction>'
write(*,*) 'was built with a different regular expression as the one defined in'
write(*,*) 'SWMM_JUNCTIONS_GUTTER_REGEX keyword. Build the file with same regular exression'
stop 'ReadDataFile - ModuleRunOff - ERR695'
endif
case default
write(*,*) 'With <BeginStormWaterGutterInteraction>/<EndStormWaterGutterInteraction>'
write(*,*) 'Use ascii file input - INITIALIZATION_METHOD keyword ASCII_FILE'
stop 'ReadDataFile - ModuleRunOff - ERR694.8'
end select
endif
call KillFillMatrix(NumberOfStormWaterNodesID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR696'
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0697'
!NumberOfStormWaterNodes points can only be <= NumberOfSewerStormWaterNodes points
call VerifyStreetGutterInteraction
else
!f not block exists <BeginStormWaterGutterInteraction>/<EndStormWaterGutterInteraction>
!do not allow to filter SWMM nodes (inside SWMM wrapper)
if (FoundSWMMRegExpression /= 0) then
write(*,*)'With trying to filter SWMM nodes with SWMM_JUNCTIONS_GUTTER_REGEX '
write(*,*)'<BeginStormWaterGutterInteraction>/<EndStormWaterGutterInteraction> must be defined in'
write(*,*) 'module Runoff.'
stop 'ReadDataFile - ModuleRunOff - ERR0696.5'
endif
!If not defined in file it will be the same as NumberOfSewerStormWaterNodes
call SetMatrixValue(Me%NumberOfStormWaterNodes, Me%Size, Me%NumberOfSewerStormWaterNodes)
endif
call RewindBuffer (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR697'
!Gets Street Gutter Length in each grid cell
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
'<BeginStreetGutterLength>', &
'<EndStreetGutterLength>', BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR699'
if (BlockFound) then
call ConstructFillMatrix ( PropertyID = StreetGutterLengthID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%StreetGutterLength, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR700'
call KillFillMatrix(StreetGutterLengthID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR710'
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0710.1'
else
write(*,*)'Missing Block <BeginStreetGutterLength> / <EndStreetGutterLength>'
stop 'ReadDataFile - ModuleRunOff - ERR711'
endif
!!Get file with 1D interactions (for External 1D Model interpolation of level and integration of computed flow)
!call GetData(MappingFileName, &
! Me%ObjEnterData,iflag, &
! SearchType = FromBlock, &
! keyword = '1D_INTERACTION_MAPPING_FILE', &
! ClientModule = 'ModuleRunoff', &
! STAT = STAT_CALL)
!if (STAT_CALL .NE. SUCCESS_) &
! stop 'ReadDataFile - ModuleRunOff - ERR800'
!
!if (iflag .EQ. 1) then
! call Read1DInteractionMapping(MappingFileName)
!endif
endif
!Get mapping to river in case of external model or DN
if (Me%StormWaterModel .or. Me%ObjDrainageNetwork /= 0) then
!Get file with 1D interactions (for External 1D Model interpolation of level and integration of computed flow)
call GetData(MappingFileName, &
Me%ObjEnterData,iflag, &
SearchType = FromFile, &
keyword = '1D_INTERACTION_MAPPING_FILE', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'ReadDataFile - ModuleRunOff - ERR800'
if (iflag .EQ. 1) then
Me%Use1D2DInteractionMapping = .true.
call Read1DInteractionMapping(MappingFileName)
endif
endif
!Write Max Flow Modulus File
call GetData(Me%Output%WriteMaxFlowModulus, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'WRITE_MAX_FLOW_FILE', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR720'
if(Me%Output%WriteMaxFlowModulus) then
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%MaxFlowModulusFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR730'
Me%Output%MaxFlowModulusFile = trim(adjustl(Me%Output%MaxFlowModulusFile))//"MaxRunOff.dat"
end if
!Write all 3 flood layers
call GetData(Me%Output%OutputFloodRisk, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'OUTPUT_FLOOD_RISK', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR731'
if (Me%Output%OutputFloodRisk) then
Me%Output%WriteMaxWaterColumn = .true.
Me%Output%WriteVelocityAtMaxWaterColumn = .true.
Me%Output%WriteMaxFloodRisk = .true.
Me%Output%WriteFloodPeriod = .true.
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%MaxWaterColumnFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0732'
Me%Output%MaxWaterColumnFile = trim(adjustl(Me%Output%MaxWaterColumnFile))//"MaxWaterColumn.dat"
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%VelocityAtMaxWaterColumnFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0733'
Me%Output%VelocityAtMaxWaterColumnFile = trim(adjustl(Me%Output%VelocityAtMaxWaterColumnFile))&
&//"VelocityAtMaxWaterColumn.dat"
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%MaxFloodRiskFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0734'
Me%Output%MaxFloodRiskFile = trim(adjustl(Me%Output%MaxFloodRiskFile))//"MaxFloodRisk.dat"
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%FloodPeriodFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0740'
Me%Output%FloodPeriodFile = trim(adjustl(Me%Output%FloodPeriodFile))//"FloodPeriod.dat"
else
!Write Max water column
call GetData(Me%Output%WriteMaxWaterColumn, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'WRITE_MAX_WATER_COLUMN', &
default = .true., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR740'
if(Me%Output%WriteMaxWaterColumn) then
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%MaxWaterColumnFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0750'
Me%Output%MaxWaterColumnFile = trim(adjustl(Me%Output%MaxWaterColumnFile))//"MaxWaterColumn.dat"
!Write velocity at maximum water column
call GetData(Me%Output%WriteVelocityAtMaxWaterColumn, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'WRITE_VELOCITY_AT_MAX_WATER_COLUMN', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR760'
if(Me%Output%WriteVelocityAtMaxWaterColumn) then
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%VelocityAtMaxWaterColumnFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0770'
Me%Output%VelocityAtMaxWaterColumnFile = trim(adjustl(Me%Output%VelocityAtMaxWaterColumnFile))&
&//"VelocityAtMaxWaterColumn.dat"
end if
endif
!Write max water column * velocity
call GetData(Me%Output%WriteMaxFloodRisk, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'WRITE_MAX_FLOOD_RISK', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR780'
if(Me%Output%WriteMaxFloodRisk) then
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%MaxFloodRiskFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0790'
Me%Output%MaxFloodRiskFile = trim(adjustl(Me%Output%MaxFloodRiskFile))//"MaxFloodRisk.dat"
endif
!Write flood period
call GetData(Me%Output%WriteFloodPeriod, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'WRITE_FLOOD_PERIOD', &
default = .false., &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR792'
if(Me%Output%WriteFloodPeriod) then
!Gets the root path from the file nomfich.dat
call ReadFileName("ROOT_SRT", Me%Output%FloodPeriodFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR0793'
Me%Output%FloodPeriodFile = trim(adjustl(Me%Output%FloodPeriodFile))//"FloodPeriod.dat"
endif
endif
!factor for velocity in flood risk
if (Me%Output%OutputFloodRisk .or. Me%Output%WriteMaxFloodRisk) then
call GetData(Me%Output%FloodRiskVelCoef, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'FLOOD_RISK_VEL_COEF', &
default = 0.5, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR795'
endif
!water column limit above which the cell is considered flooded
if (Me%Output%OutputFloodRisk .or. Me%Output%WriteFloodPeriod) then
call GetData(Me%Output%FloodWaterColumnLimit, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'FLOOD_WATER_COLUMN_LIMIT', &
default = 0.05, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR796'
endif
call ReadConvergenceParameters
call ConstructTimeSeries
call StartOutputBoxFluxes
!Closes Data File
call KillEnterData (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR800'
end subroutine ReadDataFile
!--------------------------------------------------------------------------
subroutine Read1DInteractionMapping(filename)
!Arguments-------------------------------------------------------------
character(len=PathLength) :: filename
!Local----------------------------------------------------------------
integer :: mapping1DObjEnterData, ClientNumber, STAT_CALL
integer :: iflag
type(T_NodeGridPoint), pointer :: NewNodeGridPoint, NodeGridPointFlux
type(T_BankGridPoint), pointer :: NewBankGridPoint, BankGridPointFlux, BankGridPointUp, BankGridPointDown
type(T_MarginGridPoint), pointer :: NewMarginGridPoint
logical :: BlockFound, FoundFlux, FoundUp, FoundDown
!Begin----------------------------------------------------------------
mapping1DObjEnterData = 0
call ConstructEnterData(mapping1DObjEnterData, filename, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR10'
do1: do
!Constructs Node Grid Point that will have 1D river level
call ExtractBlockFromBuffer(mapping1DObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<BeginNodeGridPoint>', &
block_end = '<EndNodeGridPoint>', &
BlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL == SUCCESS_ .and. BlockFound) then
allocate (NewNodeGridPoint)
nullify(NewNodeGridPoint%Prev,NewNodeGridPoint%Next)
call GetData(NewNodeGridPoint%ID, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'ID', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR020'
call GetData(NewNodeGridPoint%GridI, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'GRID_I', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR030'
call GetData(NewNodeGridPoint%GridJ, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'GRID_J', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR040'
call AddNodeGridPoint(NewNodeGridPoint)
else
call Block_Unlock(mapping1DObjEnterData, ClientNumber, STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR050'
exit do1
endif
enddo do1
do2: do
!Constructs Bank Grid Point that will recieve 1D river level from Node Grid Point
call ExtractBlockFromBuffer(mapping1DObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<BeginBankGridPoint>', &
block_end = '<EndBankGridPoint>', &
BlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL == SUCCESS_ .and. BlockFound) then
allocate (NewBankGridPoint)
nullify(NewBankGridPoint%Prev,NewBankGridPoint%Next)
call GetData(NewBankGridPoint%ID, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'ID', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR060'
call GetData(NewBankGridPoint%GridI, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'GRID_I', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR070'
call GetData(NewBankGridPoint%GridJ, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'GRID_J', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR080'
!link to NodeGridPoint
call GetData(NewBankGridPoint%NGPId, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'NGP_ID', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR090'
call AddBankGridPoint(NewBankGridPoint)
else
call Block_Unlock(mapping1DObjEnterData, ClientNumber, STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR095'
exit do2
endif
enddo do2
do3: do
!Constructs Margin Grid Point that will have 1D river level interpolated from Bank River Points
call ExtractBlockFromBuffer(mapping1DObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<BeginMarginGridPoint>', &
block_end = '<EndMarginGridPoint>', &
BlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL == SUCCESS_ .and. BlockFound) then
allocate (NewMarginGridPoint)
nullify(NewMarginGridPoint%Prev,NewMarginGridPoint%Next)
call GetData(NewMarginGridPoint%ID, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'ID', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR0100'
call GetData(NewMarginGridPoint%GridI, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'GRID_I', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR0110'
call GetData(NewMarginGridPoint%GridJ, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'GRID_J', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR0120'
!link to BankGridPoint Upstream
call GetData(NewMarginGridPoint%BGPUpId, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'BGP_INTERPOLATION_1_ID', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR0130'
!link to BankGridPoint Downstream
call GetData(NewMarginGridPoint%BGPDownId, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'BGP_INTERPOLATION_2_ID', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR0140'
!interpolation fraction from Upstream for each MarginGridPoint
call GetData(NewMarginGridPoint%InterpolationFraction, &
mapping1DObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'INTERPOLATION_FRACTION', &
ClientModule = 'ModuleRunoff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR0150'
if (NewMarginGridPoint%InterpolationFraction .lt. 0.0 .or. NewMarginGridPoint%InterpolationFraction .gt. 1.0) then
write (*,*)
write (*,*) 'Margin GridPoint INTERPOLATION_FRACTION can not be negative or higher than 1.0'
call SetError(FATAL_, KEYWORD_, "Read1DInteractionMapping - ModuleRunOff - ERR0155")
endif
!BankGripoint to send the computed flow depends on interpolation fraction
!if less then 0.5 goes to upstream BGPm otherwise goes to downstream BGP
if (NewMarginGridPoint%InterpolationFraction .le. 0.5) then
NewMarginGridPoint%BGPIntegrateFluxId = NewMarginGridPoint%BGPUpId
else
NewMarginGridPoint%BGPIntegrateFluxId = NewMarginGridPoint%BGPDownId
endif
!associate i and j from BGP (DN) or NGP (OpenMI) to mgp to avoid searching in run-time
!BGP or NGP where to associate flux
call FindBankGridPoint(NewMarginGridPoint%BGPIntegrateFluxId, BankGridPointFlux, FoundFlux)
if (.not. FoundFlux) then
write(*,*)
write(*,*)'Not found BankGridPoint to recieve flux ', NewMarginGridPoint%BGPIntegrateFluxId
call SetError(FATAL_, KEYWORD_, "Read1DInteractionMapping - ModuleRunOff - ERR0157")
else
!in case DN flux will be interacted at Bank points. stop here
if (Me%ObjDrainageNetwork /= 0) then
NewMarginGridPoint%GridIIntegrateFlux = BankGridPointFlux%GridI
NewMarginGridPoint%GridJIntegrateFlux = BankGridPointFlux%GridJ
!in case openMI need to go to NGP (one level down) to make sure the flux is placed in correct node
elseif (Me%StormWaterModel) then
call FindNodeGridPoint(BankGridPointFlux%NGPId, NodeGridPointFlux, FoundFlux)
if (.not. FoundFlux) then
write(*,*)
write(*,*)'Not found NodeGridPoint to recieve flux ', BankGridPointFlux%NGPId
call SetError(FATAL_, KEYWORD_, "Read1DInteractionMapping - ModuleRunOff - ERR0158")
else
NewMarginGridPoint%GridIIntegrateFlux = NodeGridPointFlux%GridI
NewMarginGridPoint%GridJIntegrateFlux = NodeGridPointFlux%GridJ
endif
endif
endif
call AddMarginGridPoint(NewMarginGridPoint)
else
call Block_Unlock(mapping1DObjEnterData, ClientNumber, STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'Read1DInteractionMapping - ModuleRunoff - ERR0160'
exit do3
endif
enddo do3
end subroutine Read1DInteractionMapping
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! This subroutine adds a new NodeGridPoint to the List
subroutine AddNodeGridPoint(NewNodeGridPoint)
!Arguments-------------------------------------------------------------
type(T_NodeGridPoint), pointer :: NewNodeGridPoint
!----------------------------------------------------------------------
! Add to the NodeGridPoint List a new NodeGridPoint
if (.not.associated(Me%FirstNodeGridPoint)) then
Me%NodeGridPointNumber = 1
Me%FirstNodeGridPoint => NewNodeGridPoint
Me%LastNodeGridPoint => NewNodeGridPoint
else
NewNodeGridPoint%Prev => Me%LastNodeGridPoint
Me%LastNodeGridPoint%Next => NewNodeGridPoint
Me%LastNodeGridPoint => NewNodeGridPoint
Me%NodeGridPointNumber = Me%NodeGridPointNumber + 1
end if
end subroutine AddNodeGridPoint
!--------------------------------------------------------------------------
! This subroutine adds a new BankGridPoint to the List
subroutine AddBankGridPoint(NewBankGridPoint)
!Arguments-------------------------------------------------------------
type(T_BankGridPoint), pointer :: NewBankGridPoint
!----------------------------------------------------------------------
! Add to the BankGridPoint List a new BankGridPoint
if (.not.associated(Me%FirstBankGridPoint)) then
Me%BankGridPointNumber = 1
Me%FirstBankGridPoint => NewBankGridPoint
Me%LastBankGridPoint => NewBankGridPoint
else
NewBankGridPoint%Prev => Me%LastBankGridPoint
Me%LastBankGridPoint%Next => NewBankGridPoint
Me%LastBankGridPoint => NewBankGridPoint
Me%BankGridPointNumber = Me%BankGridPointNumber + 1
end if
end subroutine AddBankGridPoint
!--------------------------------------------------------------------------
! This subroutine adds a new MarginGridPoint to the List
subroutine AddMarginGridPoint(NewMarginGridPoint)
!Arguments-------------------------------------------------------------
type(T_MarginGridPoint), pointer :: NewMarginGridPoint
!----------------------------------------------------------------------
! Add to the MarginGridPoint List a new MarginGridPoint
if (.not.associated(Me%FirstMarginGridPoint)) then
Me%MarginGridPointNumber = 1
Me%FirstMarginGridPoint => NewMarginGridPoint
Me%LastMarginGridPoint => NewMarginGridPoint
else
NewMarginGridPoint%Prev => Me%LastMarginGridPoint
Me%LastMarginGridPoint%Next => NewMarginGridPoint
Me%LastMarginGridPoint => NewMarginGridPoint
Me%MarginGridPointNumber = Me%MarginGridPointNumber + 1
end if
end subroutine AddMarginGridPoint
!--------------------------------------------------------------------------
subroutine VerifyStreetGutterInteraction
!Arguments-------------------------------------------------------------
integer :: CHUNK, i, j
!Begin-----------------------------------------------------------------
CHUNK = ChunkJ !CHUNK_J(Me%WorkSize%JLB, Me%WorkSize%JUB)
!Dont openmp since write will be messy and this is 2D field and this is constructor
!!$OMP PARALLEL PRIVATE(I,J)
!!$OMP DO SCHEDULE(DYNAMIC, CHUNK)
do1: do j = Me%WorkSize%JLB, Me%WorkSize%JUB
do2: do i = Me%WorkSize%ILB, Me%WorkSize%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
!Check if street gutter interaction is consistent with storm water interaction
!there cant be more gutter interaction points than all interaction points
!since gutter interaction is only going to be used to drive gutter flow to the cells
!that have eligible points (e.g. only discharge gutter in pluvial nodes)
if (Me%NumberOfStormWaterNodes(i,j) > Me%NumberOfSewerStormWaterNodes(i,j)) then
write(*,*)
write(*,*) 'Error: Number Of Storm Water Nodes is higher than'
write(*,*) 'Number Of Sewer Storm WaterNodes in cell: ', i, j
stop 'VerifyStreetGutterInteraction - Module Runoff - ERR01'
endif
endif
enddo do2
enddo do1
!!$OMP END DO NOWAIT
!!$OMP END PARALLEL
end subroutine VerifyStreetGutterInteraction
!--------------------------------------------------------------------------
subroutine StartOutputBoxFluxes
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: STAT_CALL
integer :: iflag
logical :: exist, opened
character(len=StringLength), dimension(:), pointer :: FluxesOutputList
character(len=StringLength), dimension(:), pointer :: ScalarOutputList
!Begin-----------------------------------------------------------------
! This keyword have two functions if exist fluxes between boxes are compute
! and the value read is the name file where the boxes are defined
call GetData(Me%Files%BoxesFile, &
Me%ObjEnterData, iflag, &
Keyword = 'BOXFLUXES', &
SearchType = FromFile, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'Subroutine StartOutputBoxFluxes - ModuleRunoff. ERR02.'
cd6 : if (iflag .EQ. 1) then
Me%Output%BoxFluxes = .true.
inquire(FILE = Me%Files%BoxesFile, EXIST = exist)
cd4 : if (exist) then
inquire(FILE = Me%Files%BoxesFile, OPENED = opened)
cd5 : if (opened) then
write(*,* )
write(*,'(A)') 'BoxFluxesFileName = ', Me%Files%BoxesFile
write(*,* ) 'Already opened.'
stop 'Subroutine StartOutputBoxFluxes; ModuleRunoff. ERR04'
end if cd5
allocate(FluxesOutputList(1), ScalarOutputList(1))
FluxesOutputList = 'runoff_water'
ScalarOutputList = 'runoff_water'
call StartBoxDif(BoxDifID = Me%ObjBoxDif, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
BoxesFilePath = Me%Files%BoxesFile, &
FluxesOutputList = FluxesOutputList, &
ScalarOutputList = ScalarOutputList, &
WaterPoints2D = Me%ExtVar%BasinPoints, &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'Subroutine StartOutputBoxFluxes - ModuleRunoff. ERR15.'
deallocate(FluxesOutputList, ScalarOutputList)
nullify (FluxesOutputList, ScalarOutputList)
else
write(*,*)
write(*,*) 'Error dont have the file box.'
write(*,'(A)') 'BoxFileName = ', Me%Files%BoxesFile
stop 'Subroutine StartOutputBoxFluxes; ModuleRunoff. ERR03'
end if cd4
else
Me%Output%BoxFluxes = .false.
end if cd6
end subroutine StartOutputBoxFluxes
!--------------------------------------------------------------------------
subroutine ReadConvergenceParameters
!Local-----------------------------------------------------------------
integer :: STAT_CALL, &
iflag, &
MIN_WATER_COLUMN_STABILIZE_flag
real :: dummy_real
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!Find deprecated keywords in data file
!----------------------------------------------------------------------
call GetData(dummy_real, &
Me%ObjEnterData, MIN_WATER_COLUMN_STABILIZE_flag, &
SearchType = FromFile, &
keyword ='MIN_WATER_COLUMN_STABILIZE', &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR010")
if (MIN_WATER_COLUMN_STABILIZE_flag > 0) then
write (*,*) '======================================================================='
write (*,*) 'The following deprecated keywords were found in RunOff data file:'
write (*,*) ''
if (MIN_WATER_COLUMN_STABILIZE_flag > 0) &
write(*,*) 'MIN_WATER_COLUMN_STABILIZE: Use STABILIZE_MIN_WATER_COLUMN instead.'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR070")
endif
!----------------------------------------------------------------------
!Read convergence options
!----------------------------------------------------------------------
call GetData(Me%CV%Stabilize, &
Me%ObjEnterData, iflag, &
keyword = 'STABILIZE', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR080")
if (iflag <= 0) then
write(*,*) 'WARNING: Missing STABILIZE keyword in RunOff input data file.'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR081")
endif
if (Me%CV%Stabilize) then
!Maximun change of water content (in %) allowed in one time step.
call GetData(Me%CV%StabilizeFactor, &
Me%ObjEnterData, iflag, &
keyword = 'STABILIZE_FACTOR', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = 0.1, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR082")
if (Me%CV%StabilizeFactor < 0.0 .or. Me%CV%StabilizeFactor > 1.0) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR083")
call GetData(Me%CV%MinimumValueToStabilize, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'STABILIZE_MIN_WATER_COLUMN', &
default = Me%MinimumWaterColumn, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR084")
if (Me%CV%MinimumValueToStabilize < Me%MinimumWaterColumn) then
write (*,*)'Invalid Minimun Water Column to Stabilize value [STABILIZE_MIN_WATER_COLUMN]'
write (*,*)'Value must be greater than MIN_WATER_COLUMN'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR085")
endif
call GetData(dummy_real, &
Me%ObjEnterData, iflag, &
keyword = 'STABILIZE_RESTART_FACTOR', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = 0., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR086")
if (dummy_real <= 0.) then
Me%CV%MinToRestart = 0
else
call CountDomainPoints(dummy_real)
endif
call GetData(Me%CV%CheckDecreaseOnly, &
Me%ObjEnterData, iflag, &
keyword = 'CHECK_DEC_ONLY', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR087")
endif
!Number of iterations threshold for starting to ask for a lower DT
call GetData(Me%CV%MinIterations, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='MIN_ITERATIONS', &
Default = 1, &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR090")
if (Me%CV%MinIterations < 1) then
write (*,*)'Invalid Minimun Iterations value [MIN_ITERATIONS]'
write (*,*)'Value must be greater than 0'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR091")
endif
!Number of iterations threshold that causes the model to stop
call GetData(Me%CV%MaxIterations, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='MAX_ITERATIONS', &
Default = 1024, &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR100")
if (Me%CV%MaxIterations < Me%CV%MinIterations) then
write (*,*)'Invalid Maximun Iterations value [MAX_ITERATIONS]'
write (*,*)'Value must be greater than the value of MIN_ITERATIONS'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR101")
endif
!% of the maximun iterations that causes the DT to be cut to the value of one internal time step
call GetData(dummy_real, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'DT_CUT_FACTOR', &
default = 0.1, &
ClientModule = 'ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR110")
if (dummy_real <= 0.0 .or. dummy_real > 1.0) then
write (*,*)'Invalid DT Cut Factor [DT_CUT_FACTOR]'
write (*,*)'Value must be >= 0.0 and < 1.0'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR111")
endif
Me%CV%StabilizeHardCutLimit = dummy_real * Me%CV%MaxIterations
!Internal Time Step Split
call GetData(Me%CV%DTSplitFactor, &
Me%ObjEnterData, iflag, &
keyword = 'DT_SPLIT_FACTOR', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = 2.0, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadConvergenceParameters - ModuleRunOff - ERR120'
if (Me%CV%DTSplitFactor <= 1.0) then
write (*,*)'Invalid DT Split Factor [DT_SPLIT_FACTOR]'
write (*,*)'Value must be greater then 1.0'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR121")
endif
call GetData(dummy_real, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='DT_FACTOR', &
Default = 1.25, &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR130")
if (dummy_real <= 1.0) then
write (*,*)'Invalid DT Factor [DT_FACTOR]'
write (*,*)'Value must be greater then 1.0'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR131")
endif
call GetData(Me%CV%DTFactorUp, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='DT_FACTOR_UP', &
Default = dummy_real, &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR140")
if (Me%CV%DTFactorUp <= 1.0) then
write (*,*)'Invalid DT Factor Up [DT_FACTOR_UP]'
write (*,*)'Value must be greater then 1.0'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR141")
endif
call GetData(Me%CV%DTFactorDown, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='DT_FACTOR_DOWN', &
Default = dummy_real, &
ClientModule ='ModuleRunOff', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR150")
if (Me%CV%DTFactorDown <= 1.0) then
write (*,*)'Invalid DT Factor Down [DT_FACTOR_DOWN]'
write (*,*)'Value must be greater then 1.0'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR151")
endif
call GetData(Me%CV%LimitDTCourant, &
Me%ObjEnterData, iflag, &
keyword = 'LIMIT_DT_COURANT', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR180")
if (iflag <= 0) then
write(*,*) 'WARNING: Missing LIMIT_DT_COURANT keyword in RunOff input data file.'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR181")
endif
if (Me%CV%LimitDTCourant) then
!Gets Maximum allowed Courant Number
call GetData(Me%CV%MaxCourant, &
Me%ObjEnterData, iflag, &
keyword = 'MAX_COURANT', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
Default = 1.0, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModuleRunOff - ERR181")
endif
!----------------------------------------------------------------------
end subroutine ReadConvergenceParameters
!--------------------------------------------------------------------------
subroutine CheckBoundaryCells
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: CHUNK, i, j, di, dj
real :: Sum
!Begin-----------------------------------------------------------------
CHUNK = ChunkJ !CHUNK_J(Me%WorkSize%JLB, Me%WorkSize%JUB)
!$OMP PARALLEL PRIVATE(I,J,di,dj,Sum)
!$OMP DO SCHEDULE(DYNAMIC, CHUNK)
do1: do j = Me%WorkSize%JLB, Me%WorkSize%JUB
do2: do i = Me%WorkSize%ILB, Me%WorkSize%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
!Check if near a boundary point (no diagonal)
do3: do dj = -1, 1
do4: do di = -1, 1
Sum = dj + di
if ((Me%ExtVar%BasinPoints(i+di, j+dj) == 0) .and. (Sum .eq. -1 .or. Sum .eq. 1)) then
Me%BoundaryCells(i,j) = BasinPoint
exit do3
endif
enddo do4
enddo do3
endif
enddo do2
enddo do1
!$OMP END DO NOWAIT
!$OMP END PARALLEL
end subroutine CheckBoundaryCells
!--------------------------------------------------------------------------
subroutine ReadBoundaryConditions
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: iflag, STAT_CALL
logical :: BlockFound
integer :: ClientNumber
!Begin-----------------------------------------------------------------
call GetData(Me%MaxDtmForBoundary, &
Me%ObjEnterData, iflag, &
keyword = 'MAX_DTM_FOR_BOUNDARY', &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ReadBoundaryConditions - ERR362'
if (iflag == 0) then
write(*,*)'MAX_DTM_FOR_BOUNDARY must be defined in module Runoff'
stop 'ReadBoundaryConditions - ModuleRunOff - ERR0363'
endif
call GetData(Me%BoundaryMethod, &
Me%ObjEnterData, iflag, &
keyword = 'BOUNDARY_METHOD', &
Default = ComputeFlow_, &
ClientModule = 'ModuleRunOff', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModuleRunOff - ERR360'
if (Me%BoundaryMethod /= ComputeFlow_ .and. Me%BoundaryMethod /= InstantaneousFlow_) then
write(*,*)'BOUNDARY_METHOD must be or 1 - Compute Flow or 2 - Instantaneous FlowOut'
stop 'ReadBoundaryConditions - ModuleRunOff - ERR0363.5'
endif
!it will be changed to true if the block found
Me%BoundaryImposedLevelInTime = .false.
!Search for boundary block
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadBoundaryConditions - ModuleRunoff - ERR10'
!Constructs Impermeable Fraction
call ExtractBlockFromBuffer(Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<begin_boundary>', &
block_end = '<end_boundary>', &
BlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL == SUCCESS_ .and. BlockFound) then
call ReadLevelTimeSerie()
Me%BoundaryImposedLevelInTime = .true.
else
call GetData(Me%BoundaryValue, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'BOUNDARY_VALUE', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadBoundaryConditions - ModuleRunoff - ERR120'
if (iflag == 0) then
write(*,*)'if using boundary, BOUNDARY_VALUE must be defined in ModuleRunoff'
stop 'ReadBoundaryConditions - ModulePorousMedia - ERR0110'
endif
endif
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadBoundaryConditions - ModuleRunoff - ERR0130'
end subroutine ReadBoundaryConditions
!--------------------------------------------------------------------------
subroutine ReadLevelTimeSerie()
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: iflag, STAT_CALL
!Begin-----------------------------------------------------------------
call GetData(Me%ImposedLevelTS%TimeSerie%FileName, &
Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'FILENAME', &
ClientModule = 'FillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadPiezometerTimeSerie - ModuleRunoff - ERR01'
call GetData(Me%ImposedLevelTS%TimeSerie%DataColumn, &
Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'DATA_COLUMN', &
ClientModule = 'FillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadPiezometerTimeSerie - ModuleRunoff - ERR02'
call StartTimeSerieInput(Me%ImposedLevelTS%TimeSerie%ObjTimeSerie, &
Me%ImposedLevelTS%TimeSerie%FileName, &
Me%ObjTime, &
CheckDates = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadPiezometerTimeSerie - ModuleRunoff - ERR03'
end subroutine ReadLevelTimeSerie
!------------------------------------------------------------------
subroutine ModifyBoundaryLevel
!Local-----------------------------------------------------------------
! character(5) :: AuxChar
!Begin-----------------------------------------------------------------
call UpDateLevelValue(Me%ExtVar%Now)
!boundary values are given by the timeserie value everywhere
Me%BoundaryValue = Me%ImposedLevelTS%DefaultValue
end subroutine ModifyBoundaryLevel
!--------------------------------------------------------------------------
subroutine UpDateLevelValue(CurrentTime)
!Arguments-------------------------------------------------------------
type(T_Time) :: CurrentTime
!Local-----------------------------------------------------------------
logical :: TimeCycle
type (T_Time) :: Time1, Time2, InitialDate
real :: Value1, Value2
! real :: dt1, dt2
integer :: STAT_CALL
!Begin-----------------------------------------------------------------
call GetTimeSerieInitialData(Me%ImposedLevelTS%TimeSerie%ObjTimeSerie, InitialDate, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'UpDateLeveTimeSerielValue - ModuleRunoff - ERR01'
if (CurrentTime >= InitialDate) then
!Gets Value for current Time
call GetTimeSerieValue (Me%ImposedLevelTS%TimeSerie%ObjTimeSerie, &
CurrentTime, &
Me%ImposedLevelTS%TimeSerie%DataColumn, &
Time1, Value1, Time2, Value2, TimeCycle, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'UpDateLeveTimeSerielValue - ModuleRunoff - ERR10'
if (TimeCycle) then
Me%ImposedLevelTS%DefaultValue = Value1
!Piezometer%TimeSerieHasData = .true.
else
!Interpolates Value for current instant
call InterpolateValueInTime(CurrentTime, Time1, Value1, Time2, Value2, Me%ImposedLevelTS%DefaultValue)
endif
else
write(*,*) 'Piezometer time serie does not have data'
write(*,*) 'for the beggining of the simulation'
write(*,*) 'Piezometer name: ', Me%ImposedLevelTS%TimeSerie%FileName
stop 'UpDateLeveTimeSerielValue - ModuleRunoff - ERR20'
!Piezometer%TimeSerieHasData = .false.
endif
end subroutine UpDateLevelValue
!------------------------------------------------------------------
subroutine CountDomainPoints (percent)
!Arguments-------------------------------------------------------------
real :: percent
!Local-----------------------------------------------------------------
integer :: i, j
integer :: count
!Begin-----------------------------------------------------------------
count = 0
!Initializes Water Column
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
count = count + 1
endif
enddo
enddo
Me%CV%MinToRestart = max(int(count * percent), 0)
end subroutine CountDomainPoints
!-------------------------------------------------------------------------
subroutine InitializeVariables
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: i, j
integer :: di, dj
real :: lowestValue
!Begin-----------------------------------------------------------------
if (Me%PresentInitialWaterColumn) then
!Initializes Water Column
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
Me%myWaterLevel(i, j) = Me%InitialWaterColumn(i,j) + Me%ExtVar%Topography(i, j)
Me%MyWaterColumn(i, j) = Me%InitialWaterColumn(i,j)
Me%MyWaterColumnOld(i, j) = Me%MyWaterColumn(i,j)
Me%StabilityPoints(i, j) = 1
endif
enddo
enddo
elseif (Me%PresentInitialWaterLevel) then
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
if (Me%InitialWaterLevel(i,j) > Me%ExtVar%Topography(i, j)) then
Me%myWaterLevel(i, j) = Me%InitialWaterLevel(i,j)
Me%MyWaterColumn(i, j) = Me%InitialWaterLevel(i,j) - Me%ExtVar%Topography(i, j)
else
Me%myWaterLevel(i, j) = Me%ExtVar%Topography(i, j)
Me%MyWaterColumn(i, j) = 0.0
endif
Me%MyWaterColumnOld(i, j) = Me%MyWaterColumn(i,j)
Me%StabilityPoints(i, j) = 1
endif
enddo
enddo
endif
if (Me%Buildings) then
!Checks Building Height
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
if (Me%BuildingsHeight(i, j) .lt. 0.0) then
write(*,*)'Buildings Height must be greater then 0', i, j
stop 'InitializeVariables - ModuleRunOff - ERR01'
endif
endif
enddo
enddo
endif
!Finds lowest neighbor for from D8
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
!Finds lowest neighbour
lowestValue = Me%ExtVar%Topography(i, j)
do dj = -1, 1
do di = -1, 1
if (dj /= 0 .and. di /= 0 .and. Me%ExtVar%BasinPoints(i+di, j+dj) == BasinPoint) then
!Checks lowest neighbor
if (Me%ExtVar%Topography(i + di, j + dj) < lowestValue) then
lowestValue = Me%ExtVar%Topography(i + di, j + dj)
Me%LowestNeighborI(i, j) = i + di
Me%LowestNeighborJ(i, j) = j + dj
endif
endif
enddo
enddo
endif
enddo
enddo
!If drainage network module is associated and simple interaction, then don't apply stability
!to river points
if (Me%ObjDrainageNetwork /= 0 .and. Me%SimpleChannelInteraction) then
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%RiverPoints (i, j) == BasinPoint) then
Me%StabilityPoints(i, j) = 0
endif
enddo
enddo
endif
if (Me%RouteDFourPoints) then
!Checks if a given point is a DFourSink Point -> No point in the four direction is lower then the current point
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
if (((Me%ExtVar%BasinPoints(i+1, j) == BasinPoint .and. &
Me%ExtVar%Topography (i+1, j) >= Me%ExtVar%Topography(i, j)) .or. &
Me%ExtVar%BasinPoints(i+1, j) /= BasinPoint) .and. &
((Me%ExtVar%BasinPoints(i-1, j) == BasinPoint .and. &
Me%ExtVar%Topography (i-1, j) >= Me%ExtVar%Topography(i, j)) .or. &
Me%ExtVar%BasinPoints(i-1, j) /= BasinPoint) .and. &
((Me%ExtVar%BasinPoints(i, j+1) == BasinPoint .and. &
Me%ExtVar%Topography (i, j+1) >= Me%ExtVar%Topography(i, j)) .or. &
Me%ExtVar%BasinPoints(i, j+1) /= BasinPoint) .and. &
((Me%ExtVar%BasinPoints(i, j-1) == BasinPoint .and. &
Me%ExtVar%Topography (i, j-1) >= Me%ExtVar%Topography(i, j)) .or. &
Me%ExtVar%BasinPoints(i, j-1) /= BasinPoint)) then
if (Me%LowestNeighborI(i, j) /= i .or. Me%LowestNeighborJ(i, j) /= j) then
Me%DFourSinkPoint(i, j) = BasinPoint
!D 4 Sink Points are not points where stability criteria is verified
Me%StabilityPoints(i, j)= 0
endif
endif
endif
enddo
enddo
!If drainage network modules is associated, then don't apply D4 on drainage network point
if (Me%ObjDrainageNetwork /= 0) then
if (.not. Me%RouteDFourPointsOnDN) then
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
!Source Point is a DNet Point
if (Me%ExtVar%RiverPoints (i, j) == BasinPoint) then
Me%DFourSinkPoint(i, j) = 0
endif
enddo
enddo
endif
endif
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%DFourSinkPoint(i, j) == BasinPoint) then
if (Me%LowestNeighborI(i, j) /= null_int) then
!Neighbors of D 4 Sink Points are not points where stability criteria is verified
Me%StabilityPoints(Me%LowestNeighborI(i, j), Me%LowestNeighborJ(i, j)) = 0
endif
endif
enddo
enddo
endif
end subroutine InitializeVariables
!--------------------------------------------------------------------------
subroutine CheckRiverNetWorkConsistency
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: i, j
integer :: ILB, IUB, JLB, JUB, STAT_CALL
real , dimension(:, :), pointer :: ChannelsNodeLength
call GetChannelsNodeLength (Me%ObjDrainageNetwork, ChannelsNodeLength, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'CheckRiverNetWorkConsistency - ModuleRunOff - ERR01'
ILB = Me%WorkSize%ILB
IUB = Me%WorkSize%IUB
JLB = Me%WorkSize%JLB
JUB = Me%WorkSize%JUB
do j = JLB, JUB
do i = ILB, IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
if (Me%ExtVar%RiverPoints(i, j) == BasinPoint) then
if (ChannelsNodeLength(i, j) < 0.0) then
write(*,*)'Inconsistent River Network', i, j
stop 'CheckRiverNetWorkConsistency - ModuleRunOff - ERR02'
endif
else
if (ChannelsNodeLength(i, j) > 0.0) then
write(*,*)'Inconsistent River Network', i, j
stop 'CheckRiverNetWorkConsistency - ModuleRunOff - ERR03'
endif
endif
endif
enddo
enddo
call UnGetDrainageNetwork (Me%ObjDrainageNetwork, ChannelsNodeLength, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'CheckRiverNetWorkConsistency - ModuleRunOff - ERR04'
end subroutine CheckRiverNetWorkConsistency
!--------------------------------------------------------------------------
subroutine ConstructDischarges
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
character(len=StringLength) :: DischargeName
real :: CoordinateX, CoordinateY
logical :: CoordinatesON, IgnoreOK
integer :: Id, Jd, dn, DischargesNumber
integer :: STAT_CALL
type (T_Lines), pointer :: LineX
type (T_Polygon), pointer :: PolygonX
integer, dimension(:), pointer :: VectorI, VectorJ, VectorK
integer :: SpatialEmission, nCells
call Construct_Discharges(Me%ObjDischarges, Me%ObjTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR01'
call GetDischargesNumber(Me%ObjDischarges, DischargesNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR02'
do dn = 1, DischargesNumber
call GetDischargesGridLocalization(Me%ObjDischarges, dn, &
CoordinateX = CoordinateX, &
CoordinateY = CoordinateY, &
CoordinatesON = CoordinatesON, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR03'
call GetDischargesIDName (Me%ObjDischarges, dn, DischargeName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR03'
if (CoordinatesON) then
call GetXYCellZ(Me%ObjHorizontalGrid, CoordinateX, CoordinateY, Id, Jd, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR04'
if (Id < 0 .or. Jd < 0) then
call TryIgnoreDischarge(Me%ObjDischarges, dn, IgnoreOK, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR05'
if (IgnoreOK) then
write(*,*) 'Discharge outside the domain - ',trim(DischargeName),' - ',trim(Me%ModelName)
cycle
else
stop 'ModuleRunOff - ConstructDischarges - ERR06'
endif
endif
call CorrectsCellsDischarges(Me%ObjDischarges, dn, Id, Jd, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR07'
endif
!ATTENTION - NEED TO VERIFY IF DISCHARGES ARE COLLINEAR.
!Do not allow with properties since the flow used in PMP is not distributed by discharges
!and will be accounted with flow duplicating
call GetDischargeSpatialEmission(Me%ObjDischarges, dn, LineX, PolygonX, &
SpatialEmission, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR08'
if (SpatialEmission == DischPoint_) then
call GetDischargesGridLocalization(Me%ObjDischarges, dn, &
Igrid = Id, &
JGrid = Jd, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR09'
if (Me%ExtVar%BasinPoints(Id,Jd) /= WaterPoint) then
call TryIgnoreDischarge(Me%ObjDischarges, dn, IgnoreOK, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR10'
write(*,*) 'Discharge outside the domain I=',Id,' J=',Jd,'Model name=',trim(Me%ModelName)
if (IgnoreOK) then
write(*,*) 'Discharge in a land cell - ',trim(DischargeName),' - ',trim(Me%ModelName)
cycle
else
stop 'ModuleRunOff - ConstructDischarges - ERR11'
endif
endif
nCells = 1
allocate(VectorI(nCells), VectorJ(nCells))
VectorJ(nCells) = Jd
VectorI(nCells) = Id
else
if (SpatialEmission == DischLine_) then
call GetCellZInterceptByLine(Me%ObjHorizontalGrid, LineX, &
Me%ExtVar%BasinPoints, VectorI, VectorJ, VectorK, &
nCells, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR12'
if (nCells < 1) then
write(*,*) 'Discharge line intercept 0 cells'
stop 'ModuleRunOff - ConstructDischarges - ERR13'
endif
endif
if (SpatialEmission == DischPolygon_) then
call GetCellZInterceptByPolygon(Me%ObjHorizontalGrid, PolygonX, &
Me%ExtVar%BasinPoints, VectorI, VectorJ, VectorK, &
nCells, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR14'
if (nCells < 1) then
write(*,*) 'Discharge contains 0 center cells'
write(*,*) 'Or the polygon is to small and is best to a discharge in a point or'
write(*,*) 'the polygon not define properly'
stop 'ModuleRunOff - ConstructDischarges - ERR15'
endif
endif
endif
if (SpatialEmission /= DischPoint_) then
call SetLocationCellsZ (Me%ObjDischarges, dn, nCells, VectorI, VectorJ, VectorK, STAT= STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModuleRunOff - ConstructDischarges - ERR16'
!else
! if (DischVertical == DischBottom_ .or. DischVertical == DischSurf_) then
! call SetLayer (Me%ObjDischarges, dn, VectorK(nCells), STAT = STAT_CALL)
! if (STAT_CALL /= SUCCESS_) stop 'Construct_Sub_Modules - ModuleHydrodynamic - ERR220'
! endif
! deallocate(VectorI, VectorJ, VectorK)
endif
enddo
if (Me%OutPut%TimeSerieDischON) then
call Construct_Time_Serie_Discharge
endif
end subroutine ConstructDischarges
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine Construct_Time_Serie_Discharge
!Arguments-------------------------------------------------------------
!External--------------------------------------------------------------
character(len=StringLength), dimension(:), pointer :: PropertyList
!Local-----------------------------------------------------------------
integer :: STAT_CALL
integer :: dis, i, j
character(len=StringLength) :: Extension, DischargeName
!Begin-----------------------------------------------------------------
call GetDischargesNumber(Me%ObjDischarges, Me%OutPut%DischargesNumber, STAT = STAT_CALL)