Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
12045 lines (8943 sloc) 580 KB
!------------------------------------------------------------------------------
! IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE : Mohid Model
! PROJECT : Mohid Base 1
! MODULE : PorousMedia
! URL : http://www.mohid.com
! AFFILIATION : IST/MARETEC, Marine Modelling Group
! DATE : May 2003/2016
! REVISION : Frank Braunschweig - Complete Revision,
! DESCRIPTION : Simulates Water Flow in variable saturated soils
!
!------------------------------------------------------------------------------
! Keywords read in the Data File
!
! Keyword : Data Type Default !Comment
!
! BOTTOM_FILE : char - !Path to Bottom Topography File
! START_WITH_FIELD : logical 1 !Sets Theta initial Field Capacity
! CONTINUOUS : logical 0 !Continues from previous run
! STOP_ON_WRONG_DATE : logical 1 !Stops if previous run end is different from actual
! !Start
! OUTPUT_TIME : sec. sec. sec. - !Output Time
! SURFACE_OUTPUT_TIME : sec. sec. sec. - !Output Time of surface layer
! TIME_SERIE_LOCATION : char - !Path to File which defines Time Series
! CONTINUOUS_OUTPUT_FILE : logical 1 !Writes "famous" iter.log
! CONDUTIVITYFACE : integer 1 !Way to interpolate conducivity face
! !1 - Average, 2 - Maximum, 3 - Minimum, 4 - Weigthed
! HORIZONTAL_K_FACTOR : real 1.0 !Factor for Horizontal Conductivity = Kh / Kv
! CUT_OFF_THETA_LOW : real 1e-6 !Disables calculation when Theta is near ThetaR
! CUT_OFF_THETA_HIGH : real 1e-15 !Set Theta = ThetaS when Theta > ThetaS - CUT_OFF_THETA_HIGH
! MIN_ITER : integer 2 !Number of iterations below which the DT is increased
! MAX_ITER : integer 3 !Number of iterations above which the DT is decreased
! LIMIT_ITER : integer 50 !Number of iterations of a time step (for restart)
! THETA_TOLERANCE : real 0.001 !Converge Parameter
! INCREASE_DT : real 1.25 !Increase of DT when iter < MIN_ITER
! DECREASE_DT : real 0.70 !Decrease of DT when iter > MAX_ITER
!
!
!
!<beginproperty>
! NAME : Theta / waterlevel
!
! see Module FillMatrix for more options
!
!<endproperty>
!------------------------------------------------------------------------------
!
!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 ModulePorousMedia
use ModuleGlobalData
use ModuleStopWatch
use ModuleFunctions
use ModuleTime
use ModuleHDF5
use ModuleEnterData
use ModuleProfile, only : StartProfile, WriteProfile, KillProfile
use ModuleGridData, only : ConstructGridData, GetGridData, UngetGridData, &
KillGridData
use ModuleTimeSerie, only : StartTimeSerie, WriteTimeSerie, KillTimeSerie, &
GetNumberOfTimeSeries, GetTimeSerieLocation, &
TryIgnoreTimeSerie, CorrectsCellsTimeSerie, &
GetTimeSerieName, StartTimeSerieInput, &
GetTimeSerieInitialData, GetTimeSerieValue
use ModuleHorizontalGrid, only : GetHorizontalGrid, GetGridCellArea, &
WriteHorizontalGrid, UnGetHorizontalGrid, &
GetXYCellZ, GetCoordTypeList, GetGridCoordType, &
GetGridLatitudeLongitude, GetCellZInterceptByLine, &
GetCellZInterceptByPolygon
use ModuleBasinGeometry, only : GetBasinPoints, GetRiverPoints, GetCellSlope, &
UnGetBasin
use ModuleGeometry, only : ConstructGeometry, GetGeometrySize, &
GetGeometryDistances, GetGeometryKFloor, &
UnGetGeometry, ComputeInitialGeometry, &
ComputeVerticalGeometry, GetGeometryVolumes, &
GetGeometryAreas, KillGeometry
use ModuleMap, only : ConstructMap, GetWaterPoints3D, GetOpenPoints3D, &
GetComputeFaces3D, UnGetMap, &
UpdateComputeFaces3D, KillMap
use ModuleBoxDif, only : StartBoxDif, GetBoxes, GetNumberOfBoxes, UngetBoxDif, &
BoxDif, KillBoxDif
use ModuleFillMatrix, only : ConstructFillMatrix, ModifyFillMatrix, &
KillFillMatrix, GetIfMatrixRemainsConstant
use ModuleDrainageNetwork, only : GetChannelsWaterLevel, GetChannelsBottomLevel, &
GetChannelsBottomWidth, GetChannelsOpenProcess, &
GetChannelsNodeLength, UnGetDrainageNetwork
use ModuleTriangulation, only : InterPolation, ConstructTriangulation, &
SetHeightValues, GetNumberOfTriangles, &
GetTriangleList, GetNumberOfNodes, GetNodesList, &
KillTriangulation
use ModuleDischarges ,only : Construct_Discharges, GetDischargesNumber, &
GetDischargesGridLocalization, &
GetDischargeWaterFlow, GetDischargesIDName, &
TryIgnoreDischarge, GetDischargeSpatialEmission, &
CorrectsCellsDischarges, Kill_Discharges, &
SetLocationCellsZ, SetLayer, GetDischargeON, &
GetDischargeFlowDistribuiton, UnGetDischarges
use ModuleDrawing
implicit none
private
!Subroutines---------------------------------------------------------------
!Constructor
public :: ConstructPorousMedia
private :: AllocateInstance
private :: ReadDataFile
private :: ConstructBottomTopography
private :: AllocateVariables
private :: ReadSoilTypes
private :: InitialFields
!private :: ReadInitialFile
private :: ConstructHDF5Output
private :: ConstructTimeSerie
!Selector
public :: GetNextPorousMediaDT
public :: GetPotentialInfiltration
public :: GetInfiltration
public :: GetEfectiveEVTP
public :: GetGWFlowToChannels
public :: GetGWFlowToChannelsByLayer
public :: GetGWFlowOption
public :: GetGWToChannelsLayers
public :: GetGWLayer
public :: GetGWLayerOld
public :: GetFlowDischarge
public :: GetPMTotalDischargeFlowVolume
public :: GetPorousMediaTotalStoredVolume
public :: GetFluxU
public :: GetFluxV
public :: GetFluxW
public :: GetUnsatU
public :: GetUnsatV
public :: GetUnsatW
public :: GetUnsatWFinal
! public :: GetWaterColumn
public :: GetWaterContent
public :: GetRelativeWaterContent
public :: GetHead
public :: GetThetaR
public :: GetThetaS
public :: GetThetaF
public :: GetOldWaterContent
public :: GetThetaField
public :: GetComputeSoilField
public :: GetLimitThetaLow
public :: GetUnsatK
public :: GetEvaporation
public :: GetTranspiration
public :: GetEVTPVolumes
public :: GetPMBoundaryFlowVolume
public :: GetPMStoredVolume
public :: GetIgnoreWaterColumnOnEVAP
public :: GetBoundaryImposed
public :: GetBoundaryFluxWalls
public :: GetBoundaryFluxBottom
public :: GetBoundaryCells
public :: GetWaterContentForHead
public :: GetPMObjMap
public :: UnGetPorousMedia
!Modifier
public :: ModifyPorousMedia
private :: VariableSaturatedFlow
private :: Condutivity_Face
private :: CondutivityAverage
private :: CondutivityMaximum
private :: CondutivityMinimum
private :: CondutivityGeometricAverage
private :: SoilWaterVelocity
private :: SoilParameters
private :: VerticalContinuity
private :: CheckStability
private :: ComputeNextDT
! private :: ExchangeWithDrainageNetwork
private :: PorousMediaOutput
private :: OutPutTimeSeries
private :: CalculateTotalStoredVolume
!Destructor
public :: KillPorousMedia
private :: DeAllocateInstance
!private :: WriteFinalFile
!Management
private :: Ready
private :: LocateObjPorousMedia
!Interfaces----------------------------------------------------------------
interface UnGetPorousMedia
module procedure UnGetPorousMedia_R4
module procedure UnGetPorousMedia_R8
module procedure UnGetPorousMedia_R
module procedure UnGetPorousMedia_R1
module procedure UnGetPorousMedia_RI
module procedure UnGetPorousMedia_AI
module procedure UnGetPorousMedia_AI2D
end interface UnGetPorousMedia
!Parameter-------------------------------------------------------------------
!Conductivity Face
integer, parameter :: Average = 1
integer, parameter :: Maximum = 2
integer, parameter :: Minimum = 3
integer, parameter :: Weighted = 4
integer, parameter :: GeometricAvg = 5
!Evapotranspiration Method
integer, parameter :: SingleEvapotranspiration = 1
integer, parameter :: SeparateEvapotranspiration = 2
!Min Thickness UG Watercolumn
real, parameter :: MinUGThickness = 0.10
!Boundary Conditions interpolation methods
integer, parameter :: Triangulation_ = 1
integer, parameter :: InverseWeight_ = 2
integer, parameter :: NoInterpolation_ = 3
!Boundary Piezometers
character(len=9 ), parameter :: FromTimeSerie = 'TIMESERIE'
character(len=12), parameter :: SingleValue = 'CONSTANT'
!Water Contents
character(LEN = StringLength), parameter :: char_Theta = trim(adjustl('Theta' ))
!Conductivity
character(LEN = StringLength), parameter :: char_SoilID = trim(adjustl('SoilID'))
!Waterlevel
character(LEN = StringLength), parameter :: char_waterlevel = trim(adjustl('waterlevel' ))
!Infiltration conductivity - for tests only - by default nothing changes
integer, parameter :: SatCond_ = 1
integer, parameter :: UnSatCond_ = 2
!Bottom Boundary Condition
integer, parameter :: NullGradient_ = 1 !Boundary theta is bottom Theta (Null in Theta)
!and no hidrostatic pressure (always moving)
!so velocity is conductivity (Final Head gradient is 1)
!Drainage Network link formulations - gone to Global data, because of basin use?
!integer, parameter :: GWFlowToChanByCell_ = 1
!integer, parameter :: GWFlowToChanByLayer_ = 2
!GW Flow by Cell - Area method
integer, parameter :: GWFlowAreaWetPerimeter_ = 1
integer, parameter :: GWFlowAreaWetPerimeterAndAquifer_ = 2
!Restart fiels format
integer, parameter :: BIN_ = 1
integer, parameter :: HDF_ = 2
!Types---------------------------------------------------------------------
type T_OutPut
type (T_Time), pointer, dimension(:) :: OutTime => null()
type (T_Time), dimension(:), pointer :: RestartOutTime => null()
type (T_Time), dimension(:), pointer :: SurfaceOutTime => null()
integer :: NextOutPut = null_int
logical :: Yes = .false.
logical :: TimeSerieON = .false.
logical :: ProfileON = .false.
logical :: WriteRestartFile = .false.
logical :: SurfaceOutput = .false.
logical :: RestartOverwrite = .false.
logical :: BoxFluxes = .false.
integer :: NextRestartOutput = 1
integer :: NextSurfaceOutput = 1
integer :: RestartFormat = BIN_
end type T_OutPut
type T_IntegrationByHorizon
character(len=StringLength) :: Name = null_str
integer :: StartLayer, &
EndLayer
real, dimension(:,:), pointer :: Old2D => null(), &
Field2D => null()
real, dimension(:,:,:), pointer :: Old3D => null(), &
Field3D => null()
end type T_IntegrationByHorizon
type T_IntegrationInfo
logical :: Yes = .false., &
ByLayer = .false., &
ByHorizon = .false.
real, dimension(:,:), pointer :: Old2D => null(), &
Field2D => null()
real, dimension(:,:,:), pointer :: Old3D => null(), &
Field3D => null()
integer :: HorizonsCount = 0
type(T_IntegrationByHorizon), dimension(:), pointer :: Horizons => null()
end type T_IntegrationInfo
type T_IntegrationOutput
type (T_Time), dimension(:), pointer :: OutTime => null()
integer :: NextOutPut = null_int
logical :: Yes = .false., &
Initialize = .true.
type (T_IntegrationInfo) :: WaterContent, &
RelativeWaterContent, &
WaterTable, &
Infiltration, &
BoundaryBottom
real :: AccTime = 0.0, &
OldAccTime = 0.0
end type T_IntegrationOutput
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) :: BottomFile = null_str
character(PathLength) :: ASCFile = null_str
character(PathLength) :: BoxesFile = null_str
character(PathLength) :: IntegrationHDFFile = null_str
integer :: AsciiUnit = null_int
end type T_Files
type T_ExtVar
integer, dimension(:,:), pointer :: BasinPoints => null()
!ObjGeometry
real , pointer, dimension(:,: ) :: DUX => null()
real , pointer, dimension(:,: ) :: DVY => null()
real , pointer, dimension(:,: ) :: DZX => null()
real , pointer, dimension(:,: ) :: DZY => null()
real , pointer, dimension(:,:,:) :: DZZ => null()
real , pointer, dimension(:,:,:) :: DWZ => null()
real , pointer, dimension(:,: ) :: YY_IE => null()
real , pointer, dimension(:,: ) :: XX_IE => null()
real , pointer, dimension(:,:,:) :: CenterCell => null()
real , pointer, dimension(:,: ) :: Area => null()
real , pointer, dimension(:,: ) :: Topography => null()
real , pointer, dimension(:,: ) :: BottomTopoG => null()
real , pointer, dimension(:,:,:) :: AreaU => null()
real , pointer, dimension(:,:,:) :: AreaV => null()
integer, dimension(:,:), pointer :: RiverPoints => null()
integer, dimension(:,:), pointer :: KFloor => null()
real(8), pointer, dimension(:,:,:) :: CellVolume => null()
real , pointer, dimension(:,:,:) :: SZZ => null()
!Map
integer, pointer, dimension(:,:,:) :: WaterPoints3D => null()
integer, pointer, dimension(:,:,:) :: OpenPoints3D => null()
integer, pointer, dimension(:,:,:) :: ComputeFacesU3D => null()
integer, pointer, dimension(:,:,:) :: ComputeFacesV3D => null()
integer, pointer, dimension(:,:,:) :: ComputeFacesW3D => null()
real(8), dimension(:,: ), pointer :: InfiltrationColumn => null()
real, dimension(:,:,:), pointer :: TranspirationFlux => null()
real, dimension(:,: ), pointer :: PotentialEvaporationFlux => null()
real, dimension(:,: ), pointer :: ImposedInfiltration => null()
logical :: ConstructEvaporation = .false.
logical :: ConstructTranspiration = .false.
!Time
type (T_Time) :: Now
real :: DT = null_real
end type T_ExtVar
!type T_PointF
! real :: X = null_real
! real :: Y = null_real
!end type T_PointF
type T_FromTimeSerie
integer :: ObjTimeSerie = 0
character(len=StringLength) :: FileName = null_str
integer :: DataColumn = null_int
end type T_FromTimeSerie
type T_Piezometer
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
type(T_Piezometer), pointer :: Next
end type T_Piezometer
type T_InverseWeight
real :: MaxDistance = null_real
real :: IWDn = null_real
end type T_InverseWeight
type T_Nodes
real, dimension(:), pointer :: X
real, dimension(:), pointer :: Y
real, dimension(:), pointer :: Z
end type T_Nodes
type T_Triangulation
logical :: FillOutsidePoints = .false.
logical :: OutputTriangles = .false.
integer :: Instant = null_int
character(len=StringLength) :: FileName = null_str
type (T_Nodes) :: Nodes
end type T_Triangulation
type T_Boundary
logical :: ImposedLevelConstant = .false.
logical :: ImposedLevelInTime = .false.
integer :: InterpolationMethod = null_int
integer :: NumberOfPiezometers = null_int
integer :: ImposeBoundaryBottomCondition = null_int
real :: BoundaryValue = null_real
real :: MaxDtmForBoundary = null_real !Max DTM to apply Walls Boundary
real :: MinThetaFForBoundary = null_real !Min ThetaF to apply Bottom Boundary
type(T_PointF), dimension(:,:), pointer :: GridPoint => null()
real, dimension(:,:), pointer :: ImposedBoundaryLevel => null()
integer, dimension(:,:), pointer :: BoundaryCells => null()
type (T_Triangulation) :: Triangulation
type (T_InverseWeight) :: InverseWeight
type (T_Piezometer), pointer :: FirstPiezometer
end type T_Boundary
!Unsaturated Zone Types
type T_SoilOptions
logical :: CalcHorizontal = .false.
logical :: CalcDrainageNetworkFlux = .false.
integer :: CondutivityFace = null_int
logical :: Continuous = .false.
logical :: StopOnWrongDate = .false.
logical :: CheckGlobalMass = .false.
logical :: StartWithFieldCapacity = .false.
logical :: ComputeSoilField = .false.
real :: HCondFactor = null_real
real :: FCHCondFactor = null_real
logical :: LimitEVAPWaterVelocity = .false.
logical :: LimitEVAPHead = .false.
logical :: IgnoreWaterColumnOnEvap = .false.
real :: HeadLimit = null_real
integer :: DNLink = null_int
integer :: DNLinkAreaMethod = null_int
logical :: ComputeHydroPressure = .false.
integer :: InfiltrationConductivity = null_int
logical :: DryChannelsCompletely = .false.
logical :: ImposeBoundary = .false. !Boundary imposed
logical :: ImposeBoundaryWalls = .false. !Boundary imposed in wall
logical :: ImposeBoundaryBottom = .false. !Boundary imposed in bottom
logical :: WriteLog = .false.
logical :: Discharges = .false.
end type T_SoilOptions
type T_SoilType
real :: ThetaR = null_real
real :: ThetaS = null_real
real :: nfit = null_real
real :: mfit = null_real
real :: alfa = null_real
real :: lfit = null_real
real :: SatK = null_real
real :: OverSatSlope = null_real
end type T_SoilType
type T_Retention !Main parameters in the Mualem-van Genuchten retention and conductivity cuves
real, dimension(:,:,:), allocatable:: ThetaR !Minimum water content
real, dimension(:,:,:), allocatable:: ThetaS !Saturated water content
real, dimension(:,:,:), allocatable:: ThetaF !(Theta-ThetaR)/(ThetaS-ThetaR)
end type T_Retention
type T_Converge
integer :: MinIterations = 1
integer :: MaxIterations = 1024
logical :: Stabilize = .true.
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
real :: LimitThetaLo = 1.0e-15
real :: LimitThetaHi = 1.0e-15
real :: LimitThetaHiGWTable = 0.9995
real :: ThetaHydroCoef = 0.98
real :: VelHydroCoef = 1.0
logical :: CheckDecreaseOnly = .false.
real :: MinimumValueToStabilize = 0.1
integer :: MinToRestart = 0
real, allocatable, dimension(:,:,:) :: ThetaOld
real, allocatable, dimension(:,:,:) :: ThetaIni
real, allocatable, dimension(:,:,:) :: HeadIni
end type T_Converge
type T_PorousMedia
!Instaces of other Objects
integer :: InstanceID
character(len=StringLength) :: ModelName
integer :: ObjBasinGeometry = 0
integer :: ObjTime = 0
integer :: ObjGeometry = 0
integer :: ObjMap = 0
integer :: ObjHorizontalGrid = 0
integer :: ObjHorizontalMap = 0
integer :: ObjTopography = 0
integer :: ObjTimeSerie = 0
integer :: ObjHDF5 = 0
integer :: ObjIntegrationHDF5 = 0
integer :: ObjDrainageNetwork = 0
integer :: ObjBottomTopography = 0
integer :: ObjEnterData = 0
integer :: ObjProfile = 0
integer :: ObjTriangulation = 0
integer :: ObjDischarges = 0
integer :: ObjBoxDif = 0
type (T_PropertyID) :: ImpermeableFractionID
real, allocatable, dimension(:,:,:) :: ThetaField !!FieldCapacity [m3/m3]
type (T_OutPut) :: OutPut
type (T_IntegrationOutput) :: IntegrationOutput
type (T_ExtVar) :: ExtVar
type (T_Files) :: Files
type (T_Time) :: BeginTime
type (T_Time) :: EndTime
type (T_Boundary) :: Boundary
real(8), pointer, dimension(:,: ) :: WaterColumn => null()
real(8), pointer, dimension(:,: ) :: Infiltration => null()
real(8), pointer, dimension(:,: ) :: EfectiveEVTP => null()
real(8), pointer, dimension(:,: ) :: EfectiveEVAP => null()
!For Basin Water Balance
real(8) :: AccEvapFromSoil = 0.0 !m3
real(8) :: AccTranspiration = 0.0 !m3
real(8) :: AccBoundaryFlowVolume = 0.0 !m3
real(8) :: TotalDischargeFlowVolume
!Watertable Properties
real, dimension(:,:), pointer :: OldUGWaterLevel2D => null()
real, dimension(:,:), pointer :: UGWaterLevel2D => null()
real, dimension(:,:), pointer :: UGWaterDepth2D => null()
integer, dimension(:,:), pointer :: UGCell => null()
integer, dimension(:,:), pointer :: UGCell_Old => null()
real, dimension(:,:,:), allocatable :: iFlowBoundaryWalls
real, dimension(:,: ), allocatable :: iFlowBoundaryBottom
!Exchange with channels
real, dimension(:,:), allocatable :: lFlowToChannels
real, dimension(:,:,:), allocatable :: lFlowToChannelsLayer
real, dimension(:,:), allocatable :: iFlowToChannels
real, dimension(:,:,:), allocatable :: iFlowToChannelsLayer
integer, dimension(:,:), allocatable :: FlowToChannelsTopLayer
integer, dimension(:,:), allocatable :: FlowToChannelsBottomLayer
!Discharge
real, dimension(:,:,:), allocatable :: lFlowDischarge !Instantaneous Flow of discharges
real, dimension(:,:,:), allocatable :: iFlowDischarge !Integrated Flow of discharges
!Velocities
real, dimension(:,:,:), allocatable :: UnsatVelU
real, dimension(:,:,:), allocatable :: UnsatVelV
real, dimension(:,:,:), allocatable :: UnsatVelW
real, dimension(:,:,:), allocatable :: UnsatVelWFinal
!infiltration
real, pointer, dimension(:,:) :: InfiltrationVelocity => null()
real, pointer, dimension(:,:) :: ImpermeableFraction => null()
!Fluxes
real(8), dimension(:,:,:), allocatable :: FluxU
real(8), dimension(:,:,:), allocatable :: FluxV
real(8), dimension(:,:,:), allocatable :: FluxW
real(8), dimension(:,:,:), allocatable :: FluxWFinal !Flux Corrected with Vertical continuity
real, dimension(:,: ), pointer :: EvaporationFlux => null()
!Flow Properties
real, allocatable, dimension(:,:,:) :: Theta !water content on each cell [m3/m3]
real, allocatable, dimension(:,:,:) :: Head !Suction Head on each cell
real, allocatable, dimension(:,:,:) :: HydroPressure !Hydrostatic pressure
real, allocatable, dimension(:,:,:) :: FinalHead !Sum of Suction, Hydrostatic and Topography
!Common Properties
real, allocatable, dimension(:,:,:) :: SatK
integer, allocatable, dimension(:,:,:) :: SoilID
real, allocatable, dimension(:,:,:) :: UnSatK
real, allocatable, dimension(:,:,:) :: UnSatK_X
real, allocatable, dimension(:,:,:) :: UnSatK_Y
real, allocatable, dimension(:,:,:) :: UnSatK_Z
!Auxiliar SpeedUp Matrixes
logical, allocatable, dimension(:,:,:) :: CalculateHead
logical :: TranspirationExists = .false.
logical :: EvaporationExists = .false.
logical :: ImposedInfiltration = .false.
type (T_Retention ) :: RC !retention curve
type (T_Converge ) :: CV !Converge data
!Unsaturated Options
type (T_SoilOptions ) :: SoilOpt
real(8) :: TotalStoredVolume = 0.0
real(8) :: LossToGround = 0.0
!Grid size
type (T_Size3D) :: Size, WorkSize
type (T_Size2D) :: Size2D
integer :: DomainCellsNumber
!Soil Types
type (T_SoilType), dimension(:), pointer :: SoilTypes => null()
!Properties
!type (T_Property), pointer :: FirstProperty => null()
type(T_PorousMedia), pointer :: Next => null()
!Accumulated fluxes for Average flux computation (for Advection diffusion)
real(8), dimension(:,:,:), allocatable :: FluxUAcc ! => null()
real(8), dimension(:,:,:), allocatable :: FluxVAcc ! => null()
real(8), dimension(:,:,:), allocatable :: FluxWAcc ! => null()
real(8), dimension(:,:,:), allocatable :: FluxWAccFinal ! => null()
logical :: StopOnBathymetryChange = .true.
!integer :: ChunkK = 1, &
! ChunkJ = 1, &
! ChunkI = 1
end type T_PorousMedia
!Global Module Variables
type (T_PorousMedia), pointer :: FirstObjPorousMedia
type (T_PorousMedia), pointer :: Me
!--------------------------------------------------------------------------
contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine ConstructPorousMedia(ModelName, &
ObjPorousMediaID, &
ComputeTimeID, &
HorizontalGridID, &
HorizontalMapID, &
TopographyID, &
BasinGeometryID, &
DrainageNetworkID, &
CheckGlobalMass, &
ConstructEvaporation, &
ConstructTranspiration, &
GeometryID, &
MapID, &
DischargesID, &
StopOnBathymetryChange, &
STAT)
!Arguments---------------------------------------------------------------
character(len=*) :: ModelName
integer :: ObjPorousMediaID
integer :: ComputeTimeID
integer :: HorizontalGridID
integer :: HorizontalMapID
integer :: TopographyID
integer :: BasinGeometryID
integer :: DrainageNetworkID
logical :: CheckGlobalMass
logical :: ConstructEvaporation
logical :: ConstructTranspiration
logical :: StopOnBathymetryChange
integer, intent (OUT) :: GeometryID
integer, intent (OUT) :: DischargesID
integer, intent (OUT) :: MapID
integer, optional, intent(OUT) :: STAT
!External----------------------------------------------------------------
integer :: ready_
!Local-------------------------------------------------------------------
integer :: STAT_, STAT_CALL
integer :: DummyI
real :: DummyR
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mPorousMedia_)) then
nullify (FirstObjPorousMedia)
call RegisterModule (mPorousMedia_)
endif
call Ready(ObjPorousMediaID, ready_)
cd0 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
Me%ModelName = ModelName
!Associate External Instances
Me%ObjTime = AssociateInstance (mTIME_, ComputeTimeID )
Me%ObjTopography = AssociateInstance (mGRIDDATA_, TopographyID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)
Me%ObjHorizontalMap = AssociateInstance (mHORIZONTALMAP_, HorizontalMapID )
Me%ObjBasinGeometry = AssociateInstance (mBASINGEOMETRY_, BasinGeometryID )
if (DrainageNetworkID /= 0) &
Me%ObjDrainageNetwork = AssociateInstance (MDRAINAGENETWORK_, DrainageNetworkID )
Me%SoilOpt%CheckGlobalMass = CheckGlobalMass
Me%ExtVar%ConstructEvaporation = ConstructEvaporation
Me%ExtVar%ConstructTranspiration = ConstructTranspiration
!stop the model on bathymetry change?
Me%StopOnBathymetryChange = StopOnBathymetryChange
!Time
call GetComputeCurrentTime (Me%ObjTime, Me%ExtVar%Now, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructPorousMedia - ModulePorousMedia - ERR010'
call GetComputeTimeLimits (Me%ObjTime, BeginTime = Me%BeginTime, &
EndTime = Me%EndTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructPorousMedia - ModulePorousMedia - ERR020'
call GetComputeTimeStep (Me%ObjTime, Me%ExtVar%DT, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructPorousMedia - ModulePorousMedia - ERR030'
!Read data files
call ReadDataFile
call ConstructBottomTopography
!Build 3D domain
call VerticalDiscretization
call ReadConvergenceParameters
!After 3D build send to ModuleBasin IDs to be associated with ModulePorousMediaProperties and
!ModuleVegetation
GeometryID = Me%ObjGeometry
MapID = Me%ObjMap
!Updates Map (OpenPoints) for first output
call UpdateComputeFaces3D( Map_ID = Me%ObjMap, &
DummyR = DummyR, &
DummyI = DummyI, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructPorousMedia - ModulePorousMedia - ERR060'
call ReadLockExternalVar
call AllocateVariables
call ReadSoilTypes
!Build Initial fields
call InitialFields
!if (Me%SoilOpt%ImposeBoundary) call UpdateBoundaryConditions
!See if boundary piezometers exist
if (Me%SoilOpt%ImposeBoundary .and. Me%SoilOpt%ImposeBoundaryWalls) then
call CheckBoundaryCells
call ReadBoundaryConditions
!initial values
if (.not. Me%Boundary%ImposedLevelConstant) call ModifyBoundaryLevel
endif
if (Me%SoilOpt%Discharges) then
call ConstructDischarges
endif
if (Me%SoilOpt%Continuous) then
if (Me%OutPut%RestartFormat == BIN_) then
call ReadInitialFile_Bin
else if (Me%OutPut%RestartFormat == HDF_) then
call ReadInitialFile_Hdf
endif
endif
!Calculates initial Theta
if (.not. Me%SoilOpt%Continuous) then
call StartWithFieldCapacity
endif
!vegetation model growth needs field capacity computation
if (Me%SoilOpt%ComputeSoilField) then
call ComputeSoilFieldCapacity
endif
!Set initial time steps
Me%CV%CurrentDT = Me%ExtVar%DT
!Calculate initial heads
call SoilParameters
!Calculates Initial GW Cell
call CalculateUGWaterLevel
!Check if river bottom is below soil - warning
if (DrainageNetworkID /= 0) call CheckRiverBelowSoil
if (Me%OutPut%Yes .or. Me%Output%SurfaceOutput) then
call ConstructHDF5Output
endif
if (Me%IntegrationOutput%Yes) then
call ReadIntegrationConfiguration
call ConstructIntegrationHDF5Output
endif
call ConstructTimeSerie
call ConstructProfileOutput
if (Me%SoilOpt%WriteLog) call ConstructASCIIOutput
call StartOutputBoxFluxes
call KillEnterData (Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructPorousMedia - ModulePorousMedia - ERR02'
if (Me%SoilOpt%CheckGlobalMass) then
call CalculateTotalStoredVolume
endif
!First Output
if (Me%OutPut%Yes .or. Me%OutPut%SurfaceOutput) call PorousMediaOutput
call ReadUnLockExternalVar
!Returns ID
ObjPorousMediaID = Me%InstanceID
DischargesID = Me%ObjDischarges
STAT_ = SUCCESS_
else cd0
stop 'ModulePorousMedia - ConstructPorousMedia - ERR01'
end if cd0
if (present(STAT)) STAT = STAT_
!----------------------------------------------------------------------
end subroutine ConstructPorousMedia
!--------------------------------------------------------------------------
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 - ModulePorousMedia. 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; ModulePorousMedia. ERR04'
end if cd5
allocate(FluxesOutputList(1), ScalarOutputList(1))
FluxesOutputList = 'soil_water'
ScalarOutputList = 'soil_water'
call StartBoxDif(BoxDifID = Me%ObjBoxDif, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
BoxesFilePath = Me%Files%BoxesFile, &
FluxesOutputList = FluxesOutputList, &
ScalarOutputList = ScalarOutputList, &
WaterPoints3D = Me%ExtVar%WaterPoints3D, &
Size3D = Me%Size, &
WorkSize3D = Me%WorkSize, &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'Subroutine StartOutputBoxFluxes - ModulePorousMedia. 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; ModulePorousMedia. ERR03'
end if cd4
else
Me%Output%BoxFluxes = .false.
end if cd6
end subroutine StartOutputBoxFluxes
!--------------------------------------------------------------------------
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%Boundary%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 CheckRiverBelowSoil
!Local-----------------------------------------------------------------
real, dimension(:, :), pointer :: ChannelsBottomLevel
integer :: STAT_CALL, i, j
character (Len = 5) :: str_i, str_j
character (len = StringLength) :: StrWarning
!Begin-----------------------------------------------------------------
call GetChannelsBottomLevel (Me%ObjDrainageNetwork, ChannelsBottomLevel, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'CheckRiverBelowSoil - ModulePorousMedia - ERR02'
do1: do J = Me%WorkSize%JLB, Me%WorkSize%JUB
do2: do I = Me%WorkSize%ILB, Me%WorkSize%IUB
if (Me%ExtVar%RiverPoints(i, j) == OpenPoint) then
if (ChannelsBottomLevel(i,j) < Me%ExtVar%BottomTopoG(i, j)) then
!write (*,*)
!write (*,*) 'Bottom River section is lower than soil profile in cell', i,j
write(str_i, '(i4)') i
write(str_j, '(i4)') j
StrWarning = 'Bottom River section is lower than soil profile in cell (i, j): '// &
str_i//','//str_j
call SetError(WARNING_, INTERNAL_, StrWarning, OFF)
endif
endif
enddo do2
enddo do1
call UnGetDrainageNetwork (Me%ObjDrainageNetwork, ChannelsBottomLevel, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'CheckRiverBelowSoil - ModulePorousMedia - ERR06'
end subroutine CheckRiverBelowSoil
!--------------------------------------------------------------------------
subroutine AllocateInstance
!Local-----------------------------------------------------------------
type (T_PorousMedia), pointer :: NewObjPorousMedia
type (T_PorousMedia), pointer :: PreviousObjPorousMedia
!Allocates new instance
allocate (NewObjPorousMedia)
nullify (NewObjPorousMedia%Next)
!Insert New Instance into list and makes Current point to it
if (.not. associated(FirstObjPorousMedia)) then
FirstObjPorousMedia => NewObjPorousMedia
Me => NewObjPorousMedia
else
PreviousObjPorousMedia => FirstObjPorousMedia
Me => FirstObjPorousMedia%Next
do while (associated(Me))
PreviousObjPorousMedia => Me
Me => Me%Next
enddo
Me => NewObjPorousMedia
PreviousObjPorousMedia%Next => NewObjPorousMedia
endif
Me%InstanceID = RegisterNewInstance (mPorousMedia_)
end subroutine AllocateInstance
!--------------------------------------------------------------------------
subroutine ReadDataFile
!Local-----------------------------------------------------------------
integer :: STAT_CALL
integer :: iflag
!----------------------------------------------------------------------
!Reads the name of the data file from nomfich
call ReadFileName ('POROUS_DATA', Me%Files%DataFile, "PorousMedia Data File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR010'
!Reads the name of the transient HDF file from nomfich
call ReadFileName ('POROUS_HDF', Me%Files%TransientHDF, "PorousMedia HDF File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR020'
!Reads the name of the file where to store final data
call ReadFileName ('POROUS_FIN', Me%Files%FinalFile, "PorousMedia Final File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR030'
!Constructs the DataFile
call ConstructEnterData (Me%ObjEnterData, Me%Files%DataFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR050'
call GetData(Me%SoilOpt%WriteLog, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'WRITE_LOG', &
Default = .false., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR035'
if (Me%SoilOpt%WriteLog) then
!Reads the name of the file where to store ASCII data
call ReadFileName ('POROUS_ASC', Me%Files%ASCFile, "PorousMedia ITER SOL File", Me%EndTime, &
Extension = 'hyt', STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR040'
endif
!Botom file
call GetData(Me%Files%BottomFile, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'BOTTOM_FILE', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR060'
!General Options
call GetData(Me%SoilOpt%StartWithFieldCapacity, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'START_WITH_FIELD', &
Default = .true., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR070'
call GetData(Me%SoilOpt%ComputeSoilField, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'COMPUTE_SOIL_FIELD', &
Default = .false., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR080'
call GetData(Me%SoilOpt%LimitEVAPHead, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='LIMIT_EVAP_HEAD', &
Default = .false., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR100'
if (Me%SoilOpt%LimitEVAPHead) then
call GetData(Me%SoilOpt%HeadLimit, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='HEAD_LIMIT', &
Default = -100.0, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR110'
endif
call GetData(Me%SoilOpt%Continuous, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'CONTINUOUS', &
Default = .false., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR120'
if (Me%SoilOpt%Continuous) then
!Reads the name of the file where to read initial data
call ReadFileName ('POROUS_INI', Me%Files%InitialFile, "PorousMedia Initial File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR130'
call GetData(Me%SoilOpt%StopOnWrongDate, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'STOP_ON_WRONG_DATE', &
Default = .true., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR140'
endif
call GetData(Me%SoilOpt%IgnoreWaterColumnOnEvap, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='IGNORE_WATER_COLUMN_ON_EVAP', &
Default = .true., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR150'
!Output Options--------------------------------------------------------
!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
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR160'
!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 - ModulePorousMedia - ERR170'
call GetData(Me%OutPut%RestartFormat, &
Me%ObjEnterData, &
iflag, &
SearchType = FromFile, &
keyword = 'RESTART_FILE_FORMAT', &
Default = BIN_, &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR172'
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 - ModulePorousMedia - ERR175'
endif
call GetData(Me%OutPut%RestartOverwrite, &
Me%ObjEnterData, &
iflag, &
SearchType = FromFile, &
keyword = 'RESTART_FILE_OVERWRITE', &
Default = .true., &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR180'
!Output for surface output
call GetOutPutTime(Me%ObjEnterData, &
CurrentTime = Me%ExtVar%Now, &
EndTime = Me%EndTime, &
keyword = 'SURFACE_OUTPUT_TIME', &
SearchType = FromFile, &
OutPutsTime = Me%OutPut%SurfaceOutTime, &
OutPutsOn = Me%OutPut%SurfaceOutput, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR190'
!IN PROGRESS
!Sets Integrated Output Time
call GetOutPutTime(Me%ObjEnterData, &
CurrentTime = Me%ExtVar%Now, &
EndTime = Me%EndTime, &
keyword = 'INTEGRATION_TIME', &
SearchType = FromFile, &
OutPutsTime = Me%IntegrationOutput%OutTime, &
OutPutsOn = Me%IntegrationOutput%Yes, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR191'
Me%IntegrationOutput%NextOutput = 1
if (Me%IntegrationOutput%Yes) then
call ReadFileName('POROUS_INT_HDF', Me%Files%IntegrationHDFFile, &
Message = "Porous Media Integration HDF File", &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR192'
endif
!Checks consistency
if (Me%OutPut%Yes .and. Me%OutPut%SurfaceOutput) then
write(*,*)'Only normal output or 2D output can be active'
write(*,*)'OUTPUT_TIME or SURFACE_OUTPUT_TIME'
stop 'ReadDataFile - ModulePorousMedia - ERR200'
endif
!Directional Options---------------------------------------------------
call GetData(Me%SoilOpt%CalcHorizontal, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='CALC_HORIZONTAL', &
Default =.TRUE., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR210'
call GetData(Me%SoilOpt%CalcDrainageNetworkFlux, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='CALC_DRAINAGE_FLUX', &
Default =.TRUE., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR220'
! 1 - AVERAGE of the conductivity in the cells
! 2 - MAXIMUM of the conductivity in the cells
! 3 - MINIMUM of the conductivity in the cells
! 4 - WEIGTHED of the conductivity in the cells
! 5 - GEOMETRIC AVERAGE of the conductivity in the cells
call GetData(Me%SoilOpt%CondutivityFace, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='CONDUTIVITYFACE', &
Default = 1 , &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR230'
call GetData(Me%SoilOpt%HCondFactor, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='HORIZONTAL_K_FACTOR', &
Default = 1.0, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR240'
call GetData(Me%SoilOpt%FCHCondFactor, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='FC_K_FACTOR', &
Default = Me%SoilOpt%HCondFactor, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR250'
call GetData(Me%SoilOpt%LimitEVAPWaterVelocity, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'LIMIT_EVAP_WATER_VEL', &
Default = .false., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR260'
call GetData(Me%SoilOpt%DNLink, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'DN_LINK', &
Default = GWFlowToChanByCell_, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR270'
if ((Me%SoilOpt%DNLink /= GWFlowToChanByLayer_) .and. (Me%SoilOpt%DNLink /= GWFlowToChanByCell_)) then
write(*,*)' DN_LINK uncorrectly defined - 2 for GW flow to drainage network by layers '
write(*,*)' and 1 for GW flow for each cell'
stop 'ReadDataFile - ModulePorousMedia - ERR271'
endif
!New method
if (Me%SoilOpt%DNLink == GWFlowToChanByCell_) then
call GetData(Me%SoilOpt%DNLinkAreaMethod, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'DN_LINK_AREA_METHOD', &
Default = GWFlowAreaWetPerimeter_, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR275'
if ((Me%SoilOpt%DNLinkAreaMethod /= GWFlowAreaWetPerimeter_) .and. &
(Me%SoilOpt%DNLinkAreaMethod /= GWFlowAreaWetPerimeterAndAquifer_)) then
write(*,*)' DN_LINK_AREA_METHOD uncorrectly defined '
write(*,*)' 1 - area for Flux is around wet perimeter '
write(*,*)' 2 - area for Flux is around wet perimeter + aquifer if aquifer higher then river '
stop 'ReadDataFile - ModulePorousMedia - ERR276'
endif
endif
call GetData(Me%SoilOpt%ComputeHydroPressure, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'COMPUTE_HYDRO_PRESSURE', &
Default = .true., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR280'
!Impose boundary saturated levels in walls
call GetData(Me%SoilOpt%ImposeBoundaryWalls, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'IMPOSE_BOUNDARY_VALUE', &
Default = .false., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR290'
if (Me%SoilOpt%ImposeBoundaryWalls) then
Me%SoilOpt%ImposeBoundary = .true.
call GetData(Me%Boundary%MaxDtmForBoundary, &
Me%ObjEnterData, iflag, &
keyword = 'MAX_DTM_FOR_BOUNDARY', &
ClientModule = 'ModulePorousMedia', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR291'
if (iflag == 0) then
write(*,*)'MAX_DTM_FOR_BOUNDARY must be defined in module PorousMedia'
stop 'ReadDataFile - ModulePorousMedia - ERR0292'
endif
!Moved to after checking if piezometers exist
! call GetData(Me%Boundary%BoundaryValue, &
! Me%ObjEnterData, iflag, &
! SearchType = FromFile, &
! keyword = 'BOUNDARY_VALUE', &
! ClientModule ='ModulePorousMedia', &
! STAT = STAT_CALL)
! if (STAT_CALL /= SUCCESS_) stop 'GetUnSaturatedOptions - ModulePorousMedia - ERR271'
!
! if (iflag == 0) then
! write(*,*)'BOUNDARY_VALUE must be defined in module PorousMedia'
! stop 'ReadDataFile - ModulePorousMedia - ERR0230'
! endif
endif
!Impose boundary in bottom
call GetData(Me%SoilOpt%ImposeBoundaryBottom, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'IMPOSE_BOUNDARY_BOTTOM', &
Default = .false., &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR293'
if (Me%SoilOpt%ImposeBoundaryBottom) then
Me%SoilOpt%ImposeBoundary = .true.
call GetData(Me%Boundary%MinThetaFForBoundary, &
Me%ObjEnterData, iflag, &
keyword = 'MIN_THETAF_FOR_BOUNDARY', &
Default = 0.0, &
ClientModule = 'ModulePorousMedia', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR294'
call GetData(Me%Boundary%ImposeBoundaryBottomCondition, &
Me%ObjEnterData, iflag, &
keyword = 'IMPOSE_BOUNDARY_BOTTOM_CONDITION', &
Default = NullGradient_, &
ClientModule = 'ModulePorousMedia', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR295'
if (Me%Boundary%ImposeBoundaryBottomCondition /= NullGradient_) then
write(*,*)'IMPOSE_BOUNDARY_BOTTOM_CONDITION for now can only be 1 (NullGradient)'
stop 'ReadDataFile - ModulePorousMedia - ERR0296'
endif
endif
!Discharges
call GetData(Me%SoilOpt%Discharges, &
Me%ObjEnterData, iflag, &
keyword = 'DISCHARGES', &
ClientModule = 'ModulePorousMedia', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR300'
!Conductivity used for infiltration - for tests only - by default nothing changes
call GetData(Me%SoilOpt%InfiltrationConductivity, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'INFIL_CONDUCTIVITY', &
Default = SatCond_, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadDataFile - ModulePorousMedia - ERR310'
if ((Me%SoilOpt%InfiltrationConductivity .ne. SatCond_) .and. &
(Me%SoilOpt%InfiltrationConductivity .ne. UnSatCond_)) then
write(*,*)
write(*,*)'Method not known for infiltration conductivity'
write(*,*)'Please check INFIL_CONDUCTIVITY keyword'
stop 'ReadDataFile - ModulePorousMedia - ERR320'
endif
end subroutine ReadDataFile
!--------------------------------------------------------------------------
subroutine ReadConvergenceParameters
!Local-----------------------------------------------------------------
integer :: STAT_CALL, &
iflag, &
MIN_ITER_flag, &
MAX_ITER_flag, &
LIMIT_ITER_flag, &
THETA_TOLERANCE_flag, &
INCREASE_DT_flag, &
DECREASE_DT_flag, &
CUT_OFF_THETA_GW_flag
real :: dummy_real
integer :: dummy_int
!integer :: dummy_int, &
! ChunkKFactor, &
! ChunkJFactor, &
! ChunkIFactor
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!Find deprecated keywords in data file
!----------------------------------------------------------------------
call GetData(dummy_int, &
Me%ObjEnterData, MIN_ITER_flag, &
SearchType = FromFile, &
keyword ='MIN_ITER', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR010")
call GetData(dummy_int, &
Me%ObjEnterData, MAX_ITER_flag, &
SearchType = FromFile, &
keyword ='MAX_ITER', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR020")
call GetData(dummy_int, &
Me%ObjEnterData, LIMIT_ITER_flag, &
SearchType = FromFile, &
keyword ='LIMIT_ITER', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR030")
call GetData(dummy_real, &
Me%ObjEnterData, THETA_TOLERANCE_flag, &
SearchType = FromFile, &
keyword ='THETA_TOLERANCE', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR040")
call GetData(dummy_real, &
Me%ObjEnterData, INCREASE_DT_flag, &
SearchType = FromFile, &
keyword ='INCREASE_DT', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR050")
call GetData(dummy_real, &
Me%ObjEnterData, DECREASE_DT_flag, &
SearchType = FromFile, &
keyword ='DECREASE_DT', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR060")
call GetData(dummy_real, &
Me%ObjEnterData, CUT_OFF_THETA_GW_flag, &
SearchType = FromFile, &
keyword ='CUT_OFF_THETA_HIGH_GW_TABLE', &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR061")
if (MIN_ITER_flag > 0 .or. MAX_ITER_flag > 0 .or. LIMIT_ITER_flag > 0 .or. &
THETA_TOLERANCE_flag > 0 .or. INCREASE_DT_flag > 0 .or. DECREASE_DT_flag > 0 .or. &
CUT_OFF_THETA_GW_flag > 0) then
write (*,*) '======================================================================='
write (*,*) 'The following deprecated keywords were found in Porous Media data file:'
write (*,*) ''
if (MIN_ITER_flag > 0) &
write(*,*) 'MIN_ITER : Use MIN_ITERATIONS instead.'
if (MAX_ITER_flag > 0) &
write(*,*) 'MAX_ITER'
if (LIMIT_ITER_flag > 0) &
write(*,*) 'LIMIT_ITER : Use MAX_ITERATIONS instead.'
if (THETA_TOLERANCE_flag > 0) &
write(*,*) 'THETA_TOLERANCE : Use STABILIZE_FACTOR instead.'
if (INCREASE_DT_flag > 0) &
write(*,*) 'INCREASE_DT : Use DT_FACTOR_UP instead.'
if (DECREASE_DT_flag > 0) &
write(*,*) 'DECREASE_DT : Use DT_FACTOR_DOWN instead.'
if (CUT_OFF_THETA_GW_flag > 0) &
write(*,*) 'CUT_OFF_THETA_HIGH_GW_TABLE : Use GW_SAT_FACTOR instead.'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR070")
endif
!----------------------------------------------------------------------
!Read convergence options
!----------------------------------------------------------------------
!call GetData(ChunkKFactor, &
! Me%ObjEnterData, iflag, &
! keyword = 'CHUNK_K_FACTOR', &
! ClientModule = 'ModulePorousMedia', &
! Default = 3, &
! SearchType = FromFile, &
! STAT = STAT_CALL)
!if (STAT_CALL /= SUCCESS_) &
! call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR071")
!
!call GetData(ChunkJFactor, &
! Me%ObjEnterData, iflag, &
! keyword = 'CHUNK_I_FACTOR', &
! ClientModule = 'ModulePorousMedia', &
! Default = 10, &
! SearchType = FromFile, &
! STAT = STAT_CALL)
!if (STAT_CALL /= SUCCESS_) &
! call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR071")
!
!call GetData(ChunkIFactor, &
! Me%ObjEnterData, iflag, &
! keyword = 'CHUNK_J_FACTOR', &
! ClientModule = 'ModulePorousMedia', &
! Default = 10, &
! SearchType = FromFile, &
! STAT = STAT_CALL)
!if (STAT_CALL /= SUCCESS_) &
! call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR071")
Me%DomainCellsNumber = CountDomainPoints(ChunkKFactor)
call GetData(Me%CV%Stabilize, &
Me%ObjEnterData, iflag, &
keyword = 'STABILIZE', &
ClientModule = 'ModulePorousMedia', &
SearchType = FromFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR080")
if (iflag <= 0) then
write(*,*) 'WARNING: Missing STABILIZE keyword in Porous Media input data file.'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 = 'ModulePorousMedia', &
SearchType = FromFile, &
Default = 0.1, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR082")
if (Me%CV%StabilizeFactor < 0.0 .or. Me%CV%StabilizeFactor > 1.0) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR083")
call GetData(Me%CV%MinimumValueToStabilize, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'STABILIZE_MIN_FACTOR', &
default = 0.0, &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR084")
if (Me%CV%MinimumValueToStabilize < 0.0) then
write (*,*)'Invalid Minimun to Stabilize value [STABILIZE_MIN_FACTOR]'
write (*,*)'Value must be >= 0'
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - ERR086")
Me%CV%MinToRestart = max(int(Me%DomainCellsNumber * dummy_real), 0)
call GetData(Me%CV%CheckDecreaseOnly, &
Me%ObjEnterData, iflag, &
keyword = 'CHECK_DEC_ONLY', &
ClientModule = 'ModulePorousMedia', &
SearchType = FromFile, &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - 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 ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - 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 = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - 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 = 'ModulePorousMedia', &
SearchType = FromFile, &
Default = 2.0, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - ERR121")
endif
call GetData(dummy_real, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='DT_FACTOR', &
Default = 1.25, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - ERR131")
endif
call GetData(Me%CV%DTFactorUp, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='DT_FACTOR_UP', &
Default = dummy_real, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - ERR141")
endif
call GetData(Me%CV%DTFactorDown, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='DT_FACTOR_DOWN', &
Default = dummy_real, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - 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 - ModulePorousMedia - ERR151")
endif
! This value says how near the ThetaR the calculation is disconected.
! Disables calculation when Theta is near ThetaR
call GetData(Me%CV%LimitThetaLo, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='CUT_OFF_THETA_LOW', &
Default = 1.0e-15, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR160")
! This value says when Theta is converted to ThetaS
! Set Theta = ThetaS
call GetData(Me%CV%LimitThetaHi, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='CUT_OFF_THETA_HIGH', &
Default = 1.0e-15, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR170")
!avoid instabilities in searching for saturated water table. using low values as 1e-15
!usually causes variations in water table depth of order of meters from iteration to another
!because of small theta variations (speially problematic in river cells)
call GetData(Me%CV%LimitThetaHiGWTable, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='GW_SAT_FACTOR', &
Default = 0.99, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR180")
! This value says from which thetaS hydrostatic pressure is to be consider
! Set Theta = ThetaS
call GetData(Me%CV%ThetaHydroCoef, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='THETA_HYDRO_COEF', &
Default = 0.98, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR190")
call GetData(Me%CV%VelHydroCoef, &
Me%ObjEnterData, iflag, &
SearchType = FromFile, &
keyword ='VEL_HYDRO_COEF', &
Default = 1.00, &
ClientModule ='ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, KEYWORD_, "ReadConvergenceParameters - ModulePorousMedia - ERR200")
!----------------------------------------------------------------------
end subroutine ReadConvergenceParameters
!--------------------------------------------------------------------------
integer function CountDomainPoints (ChunkKFactor)
!Arguments-------------------------------------------------------------
integer :: ChunkKFactor
!Local-----------------------------------------------------------------
integer :: i, j, k
integer :: count_i, count_j, count_k, STAT_CALL
logical :: j_has, i_has
!Begin-----------------------------------------------------------------
CountDomainPoints = 0
count_i = 0
count_j = 0
count_k = 0
call GetBasinPoints (Me%ObjBasinGeometry, Me%ExtVar%BasinPoints, STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'CountDomainPoints - ModulePorousMedia - ERR010'
call GetGeometryKFloor(Me%ObjGeometry, &
Z = Me%ExtVar%KFloor, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, INTERNAL_, "CountDomainPoints; ModulePorousMedia. ERR011")
!Initializes Water Column
do j = Me%Size%JLB, Me%Size%JUB
j_has = .false.
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%BasinPoints(i, j) == BasinPoint) then
j_has = .true.
i_has = .false.
do k = Me%ExtVar%KFloor(i, j), Me%WorkSize%KUB - 1
i_has = .true.
count_k = count_k + 1
enddo
if (i_has) count_i = count_i + 1
endif
enddo
if (j_has) count_j = count_j + 1
enddo
ChunkK = max((Me%Size%KUB - Me%Size%KLB) / ChunkKFactor, 1)
CountDomainPoints = count_k
call UnGetBasin (Me%ObjBasinGeometry, Me%ExtVar%BasinPoints, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'CountDomainPoints - ModulePorousMedia - ERR020'
call UnGetGeometry(Me%ObjGeometry, Me%ExtVar%KFloor, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, INTERNAL_, "CountDomainPoints; ModulePorousMedia. ERR021")
end function CountDomainPoints
!-------------------------------------------------------------------------
subroutine ReadInitialFile_Bin
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
real :: Year_File, Month_File, Day_File
real :: Hour_File, Minute_File, Second_File
integer :: InitialFile
type (T_Time) :: BeginTime, EndTimeFile, EndTime
real :: DT_error
integer :: STAT_CALL
!----------------------------------------------------------------------
call UnitsManager(InitialFile, OPEN_FILE, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadInitialFileOld - ModulePorousMedia - ERR01'
open(Unit = InitialFile, File = Me%Files%InitialFile, Form = 'UNFORMATTED', status = 'OLD', IOSTAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadInitialFileOld - ModulePorousMedia - ERR02'
!Reads Date
read(InitialFile) Year_File, Month_File, Day_File, Hour_File, Minute_File, Second_File
call SetDate(EndTimeFile, Year_File, Month_File, Day_File, Hour_File, Minute_File, Second_File)
call GetComputeTimeLimits(Me%ObjTime, BeginTime = BeginTime, EndTime = EndTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadInitialFileOld - ModulePorousMedia - ERR03'
DT_error = EndTimeFile - BeginTime
!Avoid rounding erros
!All runs are limited to second definition - David 10-2015
!if (abs(DT_error) >= 0.01) then
if (abs(DT_error) >= 1) then
write(*,*) 'The end time of the previous run is different from the start time of this run'
write(*,*) 'Date in the file'
write(*,*) Year_File, Month_File, Day_File, Hour_File, Minute_File, Second_File
write(*,*) 'DT_error', DT_error
if (Me%SoilOpt%StopOnWrongDate) stop 'ReadInitialFileOld - ModulePorousMedia - ERR04'
endif
read(InitialFile)Me%Theta
call UnitsManager(InitialFile, CLOSE_FILE, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadInitialFileOld - ModulePorousMedia - ERR05'
end subroutine ReadInitialFile_Bin
!--------------------------------------------------------------------------
subroutine ReadInitialFile_Hdf()
!External--------------------------------------------------------------
integer :: STAT_CALL
!Local-----------------------------------------------------------------
logical :: EXIST
integer :: ILB, IUB, JLB, JUB
integer :: WorkILB, WorkIUB
integer :: WorkJLB, WorkJUB
integer :: WorkKLB, WorkKUB
integer :: ObjHDF5
integer :: HDF5_READ
type (T_Time) :: BeginTime, EndTimeFile, EndTime
real, dimension(:), pointer :: TimePointer
real :: DT_error
!----------------------------------------------------------------------
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
WorkILB = Me%WorkSize%ILB
WorkIUB = Me%WorkSize%IUB
WorkJLB = Me%WorkSize%JLB
WorkJUB = Me%WorkSize%JUB
WorkKLB = Me%WorkSize%KLB
WorkKUB = Me%WorkSize%KUB
!----------------------------------------------------------------------
inquire (FILE=trim(Me%Files%InitialFile), EXIST = Exist)
cd0: if (Exist) then
!Gets File Access Code
call GetHDF5FileAccess (HDF5_READ = HDF5_READ)
ObjHDF5 = 0
!Opens HDF5 File
call ConstructHDF5 (ObjHDF5, &
trim(Me%Files%InitialFile), &
HDF5_READ, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadInitialFile - ModulePorousMedia - ERR01'
!Get Time
call HDF5SetLimits (ObjHDF5, 1, 6, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadInitialFile - ModulePorousMedia - ERR010'
allocate(TimePointer(1:6))
call HDF5ReadData (ObjHDF5, "/Time", &
"Time", &
Array1D = TimePointer, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadInitialFile - ModulePorousMedia - ERR020'
call SetDate(EndTimeFile, TimePointer(1), TimePointer(2), &
TimePointer(3), TimePointer(4), &
TimePointer(5), TimePointer(6))
call GetComputeTimeLimits(Me%ObjTime, BeginTime = BeginTime, EndTime = EndTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadInitialFile - ModulePorousMedia - ERR030'
DT_error = EndTimeFile - BeginTime
!Avoid rounding erros - Frank 08-2001
!All runs are limited to second definition - David 10-2015
!if (abs(DT_error) >= 0.01) then
if (abs(DT_error) >= 1) then
write(*,*) 'The end time of the previous run is different from the start time of this run'
write(*,*) 'Date in the file'
write(*,*) TimePointer(1), TimePointer(2), TimePointer(3), TimePointer(4), TimePointer(5), TimePointer(6)
write(*,*) 'DT_error', DT_error
if (Me%SoilOpt%StopOnWrongDate) stop 'ReadInitialFile - ModulePorousMedia - ERR040'
endif
deallocate(TimePointer)
! Reads from HDF file the Property concentration and open boundary values
call HDF5SetLimits (ObjHDF5, WorkILB, WorkIUB, &
WorkJLB, WorkJUB, &
WorkKLB, WorkKUB, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadInitialFile - ModulePorousMedia - ERR050'
call HDF5ReadData (ObjHDF5, "/Results/water content", &
"water content", &
Array3D = GetPointer(Me%Theta), &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadInitialFile - ModulePorousMedia - ERR060'
call KillHDF5 (ObjHDF5, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadInitialFile - ModulePorousMedia - ERR070'
else
write(*,*)
stop 'ReadInitialFile - ModulePorousMedia - ERR080'
end if cd0
end subroutine ReadInitialFile_Hdf
!--------------------------------------------------------------------------
subroutine ConstructBottomTopography
!Local-----------------------------------------------------------------
integer :: STAT_CALL
!Constructs GridData
call ConstructGridData (Me%ObjBottomTopography, Me%ObjHorizontalGrid, &
FileName = Me%Files%BottomFile, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructBasin - ModuleBasin - ERR05'
end subroutine ConstructBottomTopography
!--------------------------------------------------------------------------
subroutine VerticalDiscretization
!Local-----------------------------------------------------------------
integer :: STAT_CALL
!get the topography
call GetGridData (Me%ObjTopography, Me%ExtVar%Topography, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerticalDiscretization - ModulePorousMedia - ERR00'
!Constructs GridData
call ConstructGeometry (GeometryID = Me%ObjGeometry, &
GridDataID = Me%ObjBottomTopography, &
HorizontalGridID = Me%ObjHorizontalGrid, &
HorizontalMapID = Me%ObjHorizontalMap, &
ActualTime = Me%ExtVar%Now, &
SurfaceElevation = Me%ExtVar%Topography, &
BathymTopoFactor = -1.0, &
StopOnBathymetryChange = Me%StopOnBathymetryChange, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - VerticalDiscretization - ERR01'
!Map - Soil Column
call ConstructMap (Map_ID = Me%ObjMap, &
GeometryID = Me%ObjGeometry, &
HorizontalMapID = Me%ObjHorizontalMap, &
TimeID = Me%ObjTime, &
StopOnBathymetryChange = Me%StopOnBathymetryChange, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - VerticalDiscretization - ERR02'
!Get water points
call GetWaterPoints3D (Me%ObjMap, Me%ExtVar%WaterPoints3D, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - VerticalDiscretization - ERR03'
!Geometry Size
call GetGeometrySize (Me%ObjGeometry, &
Size = Me%Size, &
WorkSize = Me%WorkSize, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerticalDiscretization - ModulePorousMedia - ERR05'
Me%Size2D%ILB = Me%Size%ILB
Me%Size2D%IUB = Me%Size%IUB
Me%Size2D%JLB = Me%Size%JLB
Me%Size2D%JUB = Me%Size%JUB
!Initial geometry
call ComputeInitialGeometry(Me%ObjGeometry, &
WaterPoints3D = Me%ExtVar%WaterPoints3D, &
SurfaceElevation = Me%ExtVar%Topography, &
ContinuesCompute = .false., &
ActualTime = Me%ExtVar%Now, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerticalDiscretization - ModulePorousMedia - ERR06'
!Checks Vertical Discretization
call GetGeometryDistances(Me%ObjGeometry, &
DWZ = Me%ExtVar%DWZ, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerticalDiscretization - ModulePorousMedia. ERR07'
call UnGetGeometry( Me%ObjGeometry, Me%ExtVar%DWZ, STAT = STAT_CALL )
if (STAT_CALL /= SUCCESS_) stop 'VerticalDiscretization - ModulePorousMedia. ERR09'
call UnGetMap(Me%ObjMap, Me%ExtVar%WaterPoints3D, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
call SetError(FATAL_, INTERNAL_, "VerticalDiscretization; ModulePorousMedia. ERR10")
call UngetGridData (Me%ObjTopography, Me%ExtVar%Topography, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerticalDiscretization - ModulePorousMedia - ERR11'
end subroutine VerticalDiscretization
!--------------------------------------------------------------------------
subroutine ConstructAsciiOutPut
!Local-----------------------------------------------------------------
integer :: status
integer :: STAT_CALL
integer :: Counter
character(LEN=4) :: Number
call UnitsManager(Me%Files%AsciiUnit, OPEN_FILE, STAT = status)
if (status /= SUCCESS_) stop "ConstructAsciiOutPut - ModulePorousMedia - ERR01"
Counter = 1
do1: do
Number = ' '
write(Number, fmt='(i4)')Counter
open(UNIT = Me%Files%AsciiUnit, &
! FILE = '..'//backslash//'res'//backslash//'iter.soi_'//trim(adjustl(Number))//'.log', &
FILE = Me%Files%ASCFile, &
STATUS = "REPLACE", &
IOSTAT = STAT_CALL)
if (STAT_CALL == SUCCESS_) then
exit do1
else
Counter = Counter + 1
end if
enddo do1
write (Me%Files%AsciiUnit, FMT=*) 'YY MM DD HH MM SS Iter Time_Step '
write (Me%Files%AsciiUnit, FMT=*) ' s '
end subroutine ConstructAsciiOutPut
!--------------------------------------------------------------------------
subroutine AllocateVariables
!Local-----------------------------------------------------------------
integer :: ILB, IUB, JLB, JUB
integer :: KLB, KUB
!Bounds
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
KUB = Me%Size%KUB
KLB = Me%Size%KLB
!Water Content---------------------------------------------------------
allocate (Me%Theta (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%Head (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%HydroPressure (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%FinalHead (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%UGWaterLevel2D (ILB:IUB,JLB:JUB))
allocate(Me%UGWaterDepth2D (ILB:IUB,JLB:JUB))
allocate(Me%UGCell (ILB:IUB,JLB:JUB))
allocate(Me%UGCell_Old (ILB:IUB,JLB:JUB))
allocate(Me%WaterColumn (ILB:IUB,JLB:JUB))
allocate(Me%Infiltration (ILB:IUB,JLB:JUB))
allocate(Me%EfectiveEVTP (ILB:IUB,JLB:JUB))
allocate(Me%ImpermeableFraction (ILB:IUB,JLB:JUB))
if (Me%SoilOpt%ComputeSoilField) then
allocate (Me%ThetaField (ILB:IUB,JLB:JUB,KLB:KUB))
endif
Me%Theta = null_real
Me%Head = null_real
Me%HydroPressure = null_real
Me%FinalHead = null_real
Me%UGWaterLevel2D = null_real
Me%UGWaterDepth2D = null_real
Me%UGCell = null_int
Me%UGCell_Old = null_int
Me%WaterColumn = null_real
Me%Infiltration = null_real
Me%EfectiveEVTP = null_real
Me%ImpermeableFraction = null_real
!Conductivities--------------------------------------------------------
allocate (Me%SatK (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%UnsatK (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%UnsatK_X (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%UnsatK_Y (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%UnsatK_Z (ILB:IUB,JLB:JUB,KLB:KUB))
!SoilID
allocate (Me%SoilID (ILB:IUB,JLB:JUB,KLB:KUB))
!Speed Up Matrtix
allocate (Me%CalculateHead (ILB:IUB,JLB:JUB,KLB:KUB))
Me%SatK = null_real
Me%UnsatK = null_real
Me%UnsatK_X = null_real
Me%UnsatK_Y = null_real
Me%UnsatK_Z = null_real
Me%SoilID = null_int
Me%CalculateHead = .false.
!Retention Curve-------------------------------------------------------
allocate (Me%RC%ThetaR (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%RC%ThetaS (ILB:IUB,JLB:JUB,KLB:KUB))
allocate (Me%RC%ThetaF (ILB:IUB,JLB:JUB,KLB:KUB))
Me%RC%ThetaR = null_real
Me%RC%ThetaS = null_real
Me%RC%ThetaF = null_real
!Converge method arrays------------------------------------------------
allocate(Me%CV%HeadIni (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%CV%ThetaIni (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%CV%ThetaOld (ILB:IUB,JLB:JUB,KLB:KUB))
!Velocities------------------------------------------------------------
allocate(Me%UnsatVelU (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%UnsatVelV (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%UnsatVelW (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%UnsatVelWFinal (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxU (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxV (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxW (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxWFinal (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxUAcc (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxVAcc (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxWAcc (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FluxWAccFinal (ILB:IUB,JLB:JUB,KLB:KUB))
if (Me%ExtVar%ConstructEvaporation) then
allocate(Me%EvaporationFlux (ILB:IUB,JLB:JUB))
allocate(Me%EfectiveEVAP (ILB:IUB,JLB:JUB))
endif
allocate(Me%InfiltrationVelocity(ILB:IUB,JLB:JUB ))
Me%UnsatVelU = 0.0
Me%UnsatVelV = 0.0
Me%UnsatVelW = 0.0
Me%UnsatVelWFinal = 0.0
Me%FluxU = 0.0
Me%FluxV = 0.0
Me%FluxW = 0.0
Me%FluxWFinal = 0.0
if (Me%ExtVar%ConstructEvaporation) then
Me%EvaporationFlux = 0.0
Me%EfectiveEVAP = null_real
endif
Me%InfiltrationVelocity = null_real
!Flow to Channel
allocate(Me%iFlowToChannels (ILB:IUB,JLB:JUB))
allocate(Me%iFlowToChannelsLayer (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%lFlowToChannels (ILB:IUB,JLB:JUB))
allocate(Me%lFlowToChannelsLayer (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%FlowToChannelsTopLayer (ILB:IUB,JLB:JUB))
allocate(Me%FlowToChannelsBottomLayer (ILB:IUB,JLB:JUB))
allocate(Me%lFlowDischarge (ILB:IUB,JLB:JUB,KLB:KUB))
allocate(Me%iFlowDischarge (ILB:IUB,JLB:JUB,KLB:KUB))
Me%lFlowDischarge = 0.0 !Sets values initially to zero, so
Me%iFlowDischarge = 0.0 !model can run without Dis
if (Me%SoilOpt%ImposeBoundary) then
if (Me%SoilOpt%ImposeBoundaryWalls) then
allocate(Me%Boundary%ImposedBoundaryLevel (ILB:IUB,JLB:JUB))
Me%Boundary%ImposedBoundaryLevel = null_real
allocate(Me%Boundary%BoundaryCells (ILB:IUB,JLB:JUB))
Me%Boundary%BoundaryCells = 0
allocate(Me%iFlowBoundaryWalls (ILB:IUB,JLB:JUB,KLB:KUB))
Me%iFlowBoundaryWalls = 0.0
endif
if (Me%SoilOpt%ImposeBoundaryBottom) then
allocate(Me%iFlowBoundaryBottom (ILB:IUB,JLB:JUB))
Me%iFlowBoundaryBottom = 0.0
endif
endif
Me%iFlowToChannels = 0.0
Me%iFlowToChannelsLayer = 0.0
Me%lFlowToChannels = 0.0
Me%lFlowToChannelsLayer = 0.0
Me%FlowToChannelsBottomLayer = Me%WorkSize%KLB
Me%FlowToChannelsTopLayer = Me%WorkSize%KUB
Me%FluxUAcc = 0.0
Me%FluxVAcc = 0.0
Me%FluxWAcc = 0.0
Me%FluxWAccFinal = 0.0
end subroutine AllocateVariables
!--------------------------------------------------------------------------
subroutine ConstructDischarges
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
character(len=StringLength) :: DischargeName
real :: CoordinateX, CoordinateY
logical :: CoordinatesON, IgnoreOK
integer :: Id, Jd, Kd, dn, DischargesNumber, nc
integer :: STAT_CALL
type (T_Lines), pointer :: LineX
type (T_Polygon), pointer :: PolygonX
integer, dimension(:), pointer :: VectorI, VectorJ, VectorK
integer :: SpatialEmission, nCells, DischVertical
call Construct_Discharges(Me%ObjDischarges, Me%ObjTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR01'
call GetDischargesNumber(Me%ObjDischarges, DischargesNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR02'
do dn = 1, DischargesNumber
call GetDischargesGridLocalization(Me%ObjDischarges, dn, &
DischVertical = DischVertical, &
KGrid = Kd, &
CoordinateX = CoordinateX, &
CoordinateY = CoordinateY, &
CoordinatesON = CoordinatesON, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR03'
call GetDischargesIDName (Me%ObjDischarges, dn, DischargeName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR03'
if (CoordinatesON) then
call GetXYCellZ(Me%ObjHorizontalGrid, CoordinateX, CoordinateY, Id, Jd, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR04'
if (Id < 0 .or. Jd < 0) then
call TryIgnoreDischarge(Me%ObjDischarges, dn, IgnoreOK, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR05'
if (IgnoreOK) then
write(*,*) 'Discharge outside the domain - ',trim(DischargeName),' - ',trim(Me%ModelName)
cycle
else
stop 'ModulePorousMedia - ConstructDischarges - ERR06'
endif
endif
call CorrectsCellsDischarges(Me%ObjDischarges, dn, Id, Jd, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR07'
endif
call GetDischargeSpatialEmission(Me%ObjDischarges, dn, LineX, PolygonX, &
SpatialEmission, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR08'
!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
if (SpatialEmission == DischPoint_) then
call GetDischargesGridLocalization(Me%ObjDischarges, dn, &
DischVertical = DischVertical, &
Igrid = Id, &
JGrid = Jd, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR09'
if (Me%ExtVar%BasinPoints(Id,Jd) /= WaterPoint) then
call TryIgnoreDischarge(Me%ObjDischarges, dn, IgnoreOK, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - 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 'ModulePorousMedia - ConstructDischarges - ERR11'
endif
endif
nCells = 1
allocate(VectorI(nCells), VectorJ(nCells), VectorK(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 'ModulePorousMedia - ConstructDischarges - ERR12'
if (nCells < 1) then
write(*,*) 'Discharge line intercept 0 cells'
stop 'ModulePorousMedia - 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 'ModulePorousMedia - 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 'ModulePorousMedia - ConstructDischarges - ERR15'
endif
endif
endif
c1: select case (DischVertical)
case (DischLayer_)
VectorK(:) = Kd
case (DischDepth_)
write(*,*) "VERTICAL DISCHARGE option not active - Depth =",DischDepth_
write(*,*) 'This option is not active'
stop 'ModulePorousMedia - ConstructDischarges - ERR170'
case (DischBottom_)
call GetGeometryKFloor(Me%ObjGeometry, &
Z = Me%ExtVar%KFloor, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ModulePorousMedia - ConstructDischarges - ERR180.'
n1: do nC =1, nCells
VectorK(nC) = Me%ExtVar%Kfloor(VectorI(nC), VectorJ(nC))
enddo n1
call UnGetGeometry(Me%ObjGeometry, &
Me%ExtVar%KFloor, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ModulePorousMedia - ConstructDischarges - ERR190.'
case (DischSurf_)
VectorK(:) = Me%WorkSize%KUB
case (DischUniform_)
!do not do nothing
case default
write(*,*) "VERTICAL DISCHARGE option not known ", DischVertical
write(*,*) "The known options are : "," Bottom=",DischBottom_," Surface=",DischSurf_, &
" Layer =",DischLayer_, " Depth =",DischDepth_, &
" Uniform=",DischUniform_
stop 'ModulePorousMedia - ConstructDischarges - ERR200'
end select c1
if (SpatialEmission /= DischPoint_) then
call SetLocationCellsZ (Me%ObjDischarges, dn, nCells, VectorI, VectorJ, VectorK, STAT= STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR210'
else
if (DischVertical == DischBottom_ .or. DischVertical == DischSurf_) then
call SetLayer (Me%ObjDischarges, dn, VectorK(nCells), STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ModulePorousMedia - ConstructDischarges - ERR220'
endif
deallocate(VectorI, VectorJ, VectorK)
endif
enddo
end subroutine ConstructDischarges
!--------------------------------------------------------------------------
subroutine ReadSoilTypes
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: STAT_CALL, nSoils
logical :: SoilTypeFound
integer :: ClientNumber, iflag, SoilID
real :: thf
!Counts the number of SoilTypes
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
nSoils = 0
doS: do
!Gets soils type
call ExtractBlockFromBuffer(Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<beginsoiltype>', &
block_end = '<endsoiltype>', &
BlockFound = SoilTypeFound, &
STAT = STAT_CALL)
SF: if (STAT_CALL == SUCCESS_ .and. SoilTypeFound) then
nSoils = nSoils + 1
else
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR02'
exit doS
end if SF
enddo doS
allocate(Me%SoilTypes(nSoils))
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
doH: do
!Gets soils type
call ExtractBlockFromBuffer(Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<beginsoiltype>', &
block_end = '<endsoiltype>', &
BlockFound = SoilTypeFound, &
STAT = STAT_CALL)
HF: if (STAT_CALL == SUCCESS_ .and. SoilTypeFound) then
!Reads ID
call GetData(SoilID, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'ID', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR10'
if (iflag /= 1) then
write(*,*)'Missing ID in Soil Type definition'
stop 'ReadSoilTypes - ModulePorousMedia - ERR15'
endif
!Reads Theta R
call GetData(Me%SoilTypes(SoilID)%ThetaR, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'THETA_R', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR20'
if (iflag /= 1) then
write(*,*)'Missing THETA_R in Soil Type definition'
stop 'ReadSoilTypes - ModulePorousMedia - ERR30'
endif
!Reads Theta S
call GetData(Me%SoilTypes(SoilID)%ThetaS, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'THETA_S', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR30'
if (iflag /= 1) then
write(*,*)'Missing THETA_S in Soil Type definition'
stop 'ReadSoilTypes - ModulePorousMedia - ERR40'
endif
!Reads NFit
call GetData(Me%SoilTypes(SoilID)%NFit, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'N_FIT', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR50'
if (iflag /= 1) then
write(*,*)'Missing N_FIT in Soil Type definition'
stop 'ReadSoilTypes - ModulePorousMedia - ERR60'
endif
!Reads LFit
call GetData(Me%SoilTypes(SoilID)%LFit, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'L_FIT', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR70'
if (iflag /= 1) then
write(*,*)'Missing L_FIT in Soil Type definition'
stop 'ReadSoilTypes - ModulePorousMedia - ERR80'
endif
!Reads Alfa
call GetData(Me%SoilTypes(SoilID)%Alfa, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'ALPHA', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR90'
if (iflag /= 1) then
write(*,*)'Missing ALPHA in Soil Type definition'
stop 'ReadSoilTypes - ModulePorousMedia - ERR100'
endif
!Reads SatK
call GetData(Me%SoilTypes(SoilID)%SatK, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'SAT_K', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR110'
if (iflag /= 1) then
write(*,*)'Missing SAT_K in Soil Type definition'
stop 'ReadSoilTypes - ModulePorousMedia - ERR120'
endif
!Calculates MFIT
Me%SoilTypes(SoilID)%MFit = 1.0 - (1.0 / Me%SoilTypes(SoilID)%NFit)
!Calculates Oversat Head Slope
thf = ThetaF_ (Me%SoilTypes(SoilID)%ThetaS - Me%CV%LimitThetaHi, SoilID)
Me%SoilTypes(SoilID)%OverSatSlope = - Head_(thf, SoilID) / Me%CV%LimitThetaHi
else
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadSoilTypes - ModulePorousMedia - ERR16'
exit doH
end if HF
enddo doH
end subroutine ReadSoilTypes
!--------------------------------------------------------------------------
subroutine ReadIntegrationPropertyConfig (Info, Is3D, BlockBegin, BlockEnd, AllocateArrays, ClientNumber)
!Arguments-------------------------------------------------------------
type (T_IntegrationInfo), pointer :: Info
logical, intent(IN) :: Is3D
character(*), intent(IN) :: BlockBegin
character(*), intent(IN) :: BlockEnd
logical, intent(IN), optional :: AllocateArrays
integer, intent(IN ) :: ClientNumber
!Local-----------------------------------------------------------------
integer :: STAT_CALL, block_i
logical :: BlockFound
integer :: iflag
integer :: ILB, IUB, JLB, JUB
integer :: KLB, KUB
logical :: AllocateArrays_
!----------------------------------------------------------------------
if (present(AllocateArrays)) then
AllocateArrays_ = AllocateArrays
else
AllocateArrays_ = .true.
endif
!Bounds
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
KUB = Me%Size%KUB
KLB = Me%Size%KLB
call ExtractBlockFromBlock (Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = BlockBegin, &
block_end = BlockEnd, &
BlockInBlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR010'
if (BlockFound) then
Info%Yes = .true.
call GetData (Info%ByLayer, Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlock, &
keyword = 'BY_LAYER', &
default = .false., &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR020'
if (Info%ByLayer .and. AllocateArrays_) then
if (Is3D) then
allocate (Info%Old3d (ILB:IUB,JLB:JUB,KLB:KUB))
Info%Old3d = 0.0
allocate (Info%Field3d (ILB:IUB,JLB:JUB,KLB:KUB))
Info%Field3d = 0.0
else
allocate (Info%Old2d (ILB:IUB,JLB:JUB))
Info%Old2d = 0.0
allocate (Info%Field2d (ILB:IUB,JLB:JUB))
Info%Field2d = 0.0
endif
endif
call GetData (Info%ByHorizon, Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlock, &
keyword = 'BY_HORIZON', &
default = .false., &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR030'
if (Info%ByHorizon) then
call GetNumberOfBlocks (Me%ObjEnterData, '<begininthorizon>', '<endinthorizon>',&
FromBlockInBlock_, &
Info%HorizonsCount, &
ClientNumber, STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR040'
if (Info%HorizonsCount <= 0) stop ''
allocate (Info%Horizons (Info%HorizonsCount))
do block_i=1, Info%HorizonsCount
call ExtractBlockFromBlockFromBlock (Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<begininthorizon>', &
block_end = '<endinthorizon>', &
BlockInBlockInBlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR050'
if (BlockFound) then
call GetData (Info%Horizons(block_i)%Name, &
Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlockInBlock, &
keyword = 'NAME', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR060'
call GetData (Info%Horizons(block_i)%StartLayer, &
Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlockInBlock, &
keyword = 'START', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR070'
call GetData (Info%Horizons(block_i)%EndLayer, &
Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlockInBlock, &
keyword = 'END', &
ClientModule = 'ModulePorousMedia', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR080'
if (AllocateArrays_) then
if (Is3D) then
allocate (Info%Horizons(block_i)%Old3D (ILB:IUB,JLB:JUB,KLB:KUB))
Info%Horizons(block_i)%Old3D = 0.0
allocate (Info%Horizons(block_i)%Field3D (ILB:IUB,JLB:JUB,KLB:KUB))
Info%Horizons(block_i)%Field3D = 0.0
else
allocate (Info%Horizons(block_i)%Old2D (ILB:IUB,JLB:JUB))
Info%Horizons(block_i)%Old2D = 0.0
allocate (Info%Horizons(block_i)%Field2D (ILB:IUB,JLB:JUB))
Info%Horizons(block_i)%Field2D = 0.0
endif
endif
!
!else
! stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR090'
endif
enddo
endif
else
Info%Yes = .false.
endif
call RewindBlock(Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadIntegrationPropertyConfig - ModulePorousMedia - ERR100'
end subroutine ReadIntegrationPropertyConfig
!--------------------------------------------------------------------------
subroutine ReadIntegrationConfiguration
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: STAT_CALL
logical :: BlockFound
integer :: ClientNumber
type (T_IntegrationInfo), pointer :: aux
integer :: ILB, IUB, JLB, JUB
integer :: KLB, KUB
!----------------------------------------------------------------------
!Bounds
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
KUB = Me%Size%KUB
KLB = Me%Size%KLB
Me%IntegrationOutput%AccTime = 0.0
call RewindBuffer (Me%ObjEnterData, STAT = STAT_CALL)
call ExtractBlockFromBuffer (Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<beginintegration>', &
block_end = '<endintegration>', &
BlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop ''
if (.not. BlockFound) &
stop ''
aux => Me%IntegrationOutput%WaterContent
call ReadIntegrationPropertyConfig (aux, &
.true., &
'<beginwatercontent>', &
'<endwatercontent>', ClientNumber = ClientNumber)
aux => Me%IntegrationOutput%RelativeWaterContent
call ReadIntegrationPropertyConfig (aux, &
.true., &
'<beginrelativewatercontent>', &
'<endrelativeatercontent>', &
AllocateArrays = .false., &
ClientNumber = ClientNumber)
aux => Me%IntegrationOutput%WaterTable
call ReadIntegrationPropertyConfig (aux, &
.false., &
'<beginwatertable>', &
'<endwatertable>', ClientNumber = ClientNumber)
if (aux%yes) &
allocate (Me%OldUGWaterLevel2D (ILB:IUB,JLB:JUB))
!Not implemented yet
aux => Me%IntegrationOutput%Infiltration
call ReadIntegrationPropertyConfig (aux, &
.false., &
'<begininfiltration>', &
'<endinfiltration>', ClientNumber = ClientNumber)
aux => Me%IntegrationOutput%BoundaryBottom
call ReadIntegrationPropertyConfig (aux, &
.false., &
'<beginboundarybottom>', &
'<endboundarybottom>', ClientNumber = ClientNumber)
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadIntegrationConfiguration - ModulePorousMedia - ERR030'
end subroutine ReadIntegrationConfiguration
!--------------------------------------------------------------------------
subroutine InitialFields
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: STAT_CALL, WLNumber
integer :: KLB, KUB, ClientNumber
logical :: HorizonFound, BlockFound
logical :: AllOK
integer :: nProps, iflag
type (T_PropertyID) :: WaterLevelID
!type (T_PropertyID) :: ImpermeableFractionID
integer, allocatable, dimension(:) :: LayerControl
integer, dimension(:,:,:), pointer :: AuxPointsToFill
real, dimension(:,:,:), pointer :: AuxSoilID
integer :: i, j, k
character(LEN = StringLength) :: string
type (T_PropertyID) :: ID
allocate(LayerControl(Me%WorkSize%KLB: Me%WorkSize%KUB))
LayerControl = 0
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
doH: do
!Gets soils horizons
call ExtractBlockFromBuffer(Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<beginhorizon>', &
block_end = '<endhorizon>', &
BlockFound = HorizonFound, &
STAT = STAT_CALL)
HF: if (STAT_CALL == SUCCESS_ .and. HorizonFound) then
!Reads lower layer of horizon
call GetData(KLB, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'KLB', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR08'
if (iflag /= 1) then
write(*,*)'Missing KLB in Horizon definition'
stop 'InitialFields - ModulePorousMedia - ERR09'
endif
!Reads upper layer of horizon
call GetData(KUB, Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'KUB', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR10'
if (iflag /= 1) then
write(*,*)'Missing KUB in Horizon definition'
stop 'InitialFields - ModulePorousMedia - ERR11'
endif
do k = KLB, KUB
if (LayerControl(k) /= 0) then
write(*,*)'Inconsistent horizon definition. Layer:',k
else
LayerControl(k) = 1
endif
enddo
!Allocates Aux Matrix
allocate(AuxPointsToFill (Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB, Me%Size%KLB:Me%Size%KUB))
do i = Me%Size%ILB, Me%Size%IUB
do j = Me%Size%JLB, Me%Size%JUB
do k = Me%Size%KLB, Me%Size%KUB
if (k >= KLB .and. k <= KUB) then
AuxPointsToFill(i, j, k) = Me%ExtVar%WaterPoints3D(i, j, k)
else
AuxPointsToFill(i, j, k) = 0
endif
enddo
enddo
enddo
nProps = 0 !Control Variable
doSP: do
!Gets properties of horizon
call ExtractBlockFromBlock(Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = '<beginproperty>', &
block_end = '<endproperty>', &
BlockInBlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL .EQ. SUCCESS_ .and. BlockFound) then
!Get Name of Property
call GetData(String, Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlock, &
keyword ='NAME', &
ClientModule ='ModulePOrousMedia', &
STAT = STAT_CALL)
select case (trim(adjustl(string)))
!Water contents
case (char_Theta )
call ConstructFillMatrix(PropertyID = ID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
GeometryID = Me%ObjGeometry, &
ExtractType = FromBlockInBlock, &
PointsToFill3D = AuxPointsToFill, &
Matrix3D = GetPointer(Me%Theta), &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR97'
call KillFillMatrix (ID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR97a'
!SoilID
case (char_SoilID)
!Construct Fill Matrix cant read integer. So allocate aux matrix here.
allocate(AuxSoilID (Me%Size%ILB:Me%Size%IUB, &
Me%Size%JLB:Me%Size%JUB, &
Me%Size%KLB:Me%Size%KUB))
AuxSoilID = null_real
call ConstructFillMatrix(PropertyID = ID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
GeometryID = Me%ObjGeometry, &
ExtractType = FromBlockInBlock, &
PointsToFill3D = AuxPointsToFill, &
Matrix3D = AuxSoilID, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR98'
call KillFillMatrix (ID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR98a'
do k = Me%Size%KLB, Me%Size%KUB
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (AuxPointsToFill(i, j, k) == 1) then
Me%SoilID(i, j, k) = NINT(AuxSoilID(i, j, k))
endif
enddo
enddo
enddo
deallocate (AuxSoilID)
!Invalid
case default
write(*,*)'Invalid Property', trim(adjustl(string))
stop 'InitialFields - ModulePorousMedia - ERR99'
end select
nProps = nProps + 1
else
if (nProps /= 2) then
write(*,*)'Soil hydraulic properties incorrected defined'
stop 'InitialFields - ModulePorousMedia - ERR15'
endif
exit doSP
endif
enddo doSP
deallocate (AuxPointsToFill)
!Fills Constant Soil Variables
else
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR16'
exit doH
end if HF
enddo doH
AllOK = .true.
do k = Me%WorkSize%KLB, Me%WorkSize%KUB
if (LayerControl(k) /= 1) then
AllOK = .false.
endif
enddo
deallocate(LayerControl)
if (.not. AllOK) then
write(*,*)'Inconsistent horizon definition.'
stop 'InitialFields - ModulePorousMedia - ERR61'
endif
do k = Me%Size%KLB, Me%Size%KUB
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%Waterpoints3D(i, j, k) == 1) then
if (Me%SoilID(i, j, k) < 0) then
write(*,*)'Soils not defined for [i,j,k]', i, j, k
AllOK = .false.
endif
endif
enddo
enddo
enddo
if (.not. AllOK) then
write(*,*)'Inconsistent soil definition.'
stop 'InitialFields - ModulePorousMedia - ERR61a'
endif
!Sets Matrixes of ThetaR and ThetaS
do k = Me%Size%KLB, Me%Size%KUB
do j = Me%Size%JLB, Me%Size%JUB
do i = Me%Size%ILB, Me%Size%IUB
if (Me%ExtVar%Waterpoints3D(i, j, k) == 1) then
Me%RC%ThetaS(i, j, k) = Me%SoilTypes(Me%SoilID(i, j, k))%ThetaS
Me%RC%ThetaR(i, j, k) = Me%SoilTypes(Me%SoilID(i, j, k))%ThetaR
Me%SatK (i, j, k) = Me%SoilTypes(Me%SoilID(i, j, k))%SatK
endif
enddo
enddo
enddo
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR70'
!Constructs Water Level
call ExtractBlockFromBuffer(Me%ObjEnterData, &
ClientNumber = WLNumber, &
block_begin = '<beginwaterlevel>', &
block_end = '<endwaterlevel>', &
BlockFound = BlockFound, &
STAT = STAT_CALL)
if (STAT_CALL == SUCCESS_) then
if (.not. BlockFound) then
write(*,*)'Missing Block <beginwaterlevel> / <endwaterlevel>'
stop 'InitialFields - ModulePorousMedia - ERR80'
endif
call ConstructFillMatrix ( PropertyID = WaterLevelID, &
EnterDataID = Me%ObjEnterData, &
TimeID = Me%ObjTime, &
HorizontalGridID = Me%ObjHorizontalGrid, &
ExtractType = FromBlock, &
PointsToFill2D = Me%ExtVar%BasinPoints, &
Matrix2D = Me%UGWaterLevel2D, &
TypeZUV = TypeZ_, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR81'
call KillFillMatrix (WaterLevelID%ObjFillMatrix, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR82'
else
stop 'InitialFields - ModulePorousMedia - ERR90'
endif
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'InitialFields - ModulePorousMedia - ERR70'
!Constructs Impermeable Fraction
call ExtractBlockFromBuffer(Me%ObjEnterData, &
ClientNumber = WLNumber, &
block_begin = '<beginimpermeablefraction>', &
block_end =