Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
12577 lines (9252 sloc) 600 KB
!------------------------------------------------------------------------------
! IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE : Mohid Model
! PROJECT : Mohid Base 1
! MODULE : FillMatrix
! URL : http://www.mohid.com
! AFFILIATION : IST/MARETEC, Marine Modelling Group
! DATE : May 2003
! REVISION : Frank Braunschweig - v4.0
! DESCRIPTION : Module to Fill Matrixes
!
!------------------------------------------------------------------------------
!
!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 ModuleFillMatrix
#ifndef _NOT_IEEE_ARITHMETIC
use ieee_arithmetic
#endif
use ModuleGlobalData
use ModuleTime
use ModuleEnterData
use ModuleFunctions, only : InterpolateValueInTime, InterpolateProfile, &
SetMatrixValue, InterpolateMatrix2DInTime, &
InterpolateMatrix3DInTime, &
InterpolateAngle2DInTime, &
InterpolateAngle3DInTime, &
InterpolateLinearyMatrix2D, &
InterpolateLinearyMatrix3D, &
AngleFromFieldToGrid, &
WaveLengthHuntsApproximation
use ModuleDrawing
use ModuleBoxDif, only : StartBoxDif, GetBoxes, GetNumberOfBoxes, &
UngetBoxDif, KillBoxDif
use ModuleGridData, only : ConstructGridData, GetGridData, UnGetGridData, &
KillGridData, GetGridDataType
use ModuleHorizontalGrid, only : GetGridAngle, GetHorizontalGridSize, &
GetGridBorderLimits, GetLatitudeLongitude, &
GetDDecompOpenBorders, &
GetDDecompParameters, &
GetDDecompWorkSize2D, GetZCoordinates, &
UnGetHorizontalGrid, RotateAngleFieldToGrid, &
RotateVectorFieldToGrid, GetCheckDistortion, &
GetGridOutBorderCartLimits, &
GetGridBorderCartPolygon, &
GetHorizontalGrid, ConstructHorizontalGrid, &
KillHorizontalGrid, GetDDecompON, GetCellRotation
use ModuleTimeSerie, only : StartTimeSerieInput, GetTimeSerieValue, &
GetTimeSerieDTForNextEvent, &
GetTimeSerieTimeOfNextDataset, &
GetTimeSerieTimeOfDataset, &
GetTimeSerieDataValues,GetTimeSerieValueForIndex,&
KillTimeSerie
use ModuleGeometry, only : GetGeometryDistances, UnGetGeometry, GetGeometrySize
use ModuleHDF5, only : ConstructHDF5, HDF5ReadData, GetHDF5GroupID, &
GetHDF5FileAccess, GetHDF5GroupNumberOfItems, &
HDF5SetLimits, GetHDF5ArrayDimensions, KillHDF5, &
GetHDF5GroupExist, GetHDF5DataSetExist
use ModuleField4D, only : ConstructField4D, GetField4DNumberOfInstants, &
GetField4DInstant, ModifyField4D, &
ModifyField4DXYZ, GetField4DGeneric4DValue, &
KillField4D
use ModuleStopWatch, only : StartWatch, StopWatch
implicit none
private
!Subroutines---------------------------------------------------------------
!Constructor
public :: ConstructFillMatrix
private :: AllocateInstance
private :: ReadOptions
private :: ConstructProfileTSDefault
private :: ConstructSpaceLayers
private :: ConstructSpaceBox
private :: ConstructSpaceASCIIFile
private :: ConstructSpaceProfile
private :: ConstructAnalyticProfile
private :: ConstructSpaceTimeSerie
private :: ConstructProfileTimeSerie
private :: ConstructHDFInput
private :: HDF5TimeInstant
private :: ReadHDF5Values2D
private :: ReadHDF5Values3D
private :: ProfileTimeSerieField
!Selector
public :: GetIfMatrixRemainsConstant
public :: GetDefaultValue
public :: GetFillMatrixDTPrediction
public :: GetHDFTimeLimits
public :: GetNumberOfInstants
public :: GetTimeInstant
public :: GetNextValueForDTPred
public :: GetValuesProcessingOptions
public :: GetMultiTimeSeries
public :: UngetFillMatrix
public :: GetAngleField
public :: GetVectorialField
public :: GetAnalyticCelerityON
public :: GetAnalyticCelerity
public :: GetFilenameHDF
!Modifier
public :: ModifyFillMatrix
public :: ModifyFillMatrixVectorial
private :: ModifySpaceTimeSerie
private :: ModifyHDFInput2D
private :: ModifyHDFInput2DTime
private :: ModifyHDFInput2DGeneric4D
private :: ModifyHDFInput3D
private :: ModifyHDFInput3DTime
private :: ModifyHDFInput3DGeneric4D
private :: ModifyProfileTimeSerie
!Destructor
public :: KillFillMatrix
private :: DeAllocateInstance
!Management
private :: Ready
private :: LocateObjFillMatrix
!Data
public :: T_Station
!Interfaces----------------------------------------------------------------
interface ConstructFillMatrix
module procedure ConstructFillMatrix2D
module procedure ConstructFillMatrix3D
module procedure ConstructFillMatrix2DVectorial
module procedure ConstructFillMatrix3DVectorial
end interface ConstructFillMatrix
!interface ModifyFillMatrix
! module procedure ModifyFillMatrix
! module procedure ModifyFillMatrixVectorial
!end interface ModifyFillMatrix
interface UngetFillMatrix
module procedure UngetFillMatrix2D
module procedure UngetFillMatrix3D
end interface UngetFillMatrix
interface GetVectorialField
module procedure GetVectorialField2D
module procedure GetVectorialField3D
end interface GetVectorialField
interface GetAngleField
module procedure GetAngleField2D
module procedure GetAngleField3D
end interface GetAngleField
interface GetDefaultValue
module procedure GetDefaultValueScalar
module procedure GetDefaultValueVectorial
end interface GetDefaultValue
!Parameter-----------------------------------------------------------------
integer, parameter :: Time_4D_ = 1
integer, parameter :: Generic_4D_ = 2
integer, parameter :: Dim2D = 2
integer, parameter :: Dim3D = 3
!Initialization Methods
integer, parameter :: Constant = 1
integer, parameter :: Layers = 2
integer, parameter :: Boxes = 3
integer, parameter :: ASCIIFile = 4
integer, parameter :: Profile = 5
integer, parameter :: ReadTimeSerie = 6
integer, parameter :: ReadHDF = 7
integer, parameter :: AnalyticProfile = 8
integer, parameter :: ProfileTimeSerie = 9
integer, parameter :: Sponge = 10
integer, parameter :: MultiTimeserie = 11
integer, parameter :: AnalyticWave = 12
!Data Processing Type (used with MultiTimeserie)
integer, parameter :: Interpolate = 1 !Will interpolate the value between two times
integer, parameter :: Accumulate = 2 !Will distribute the value in time
integer, parameter :: NoProcessing = 3 !Will use original value (exact or most right one)
!Variable from file
integer, parameter :: None = 1
!Type of analytic profiles
integer, parameter :: Linear = 1
integer, parameter :: Exponential = 2
!type of values
integer, parameter :: InterpolatedValues = 1
integer, parameter :: AccumulatedValues = 2
integer, parameter :: OriginalValues = 3
!type of values
integer, parameter :: sponge_exp_ = 1
integer, parameter :: sponge_linear_ = 2
integer, parameter :: sponge_wave_stress_dump_ = 3
!wave types
integer, parameter :: SineWaveSeaLevel_ = 1
integer, parameter :: CnoidalWaveSeaLevel_ = 2
integer, parameter :: SolitartyWaveSeaLevel_ = 3
integer, parameter :: SineWaveVelX_ = 4
integer, parameter :: CnoidalWaveVelX_ = 5
integer, parameter :: SolitartyWaveVelX_ = 6
integer, parameter :: SineWaveVelY_ = 7
integer, parameter :: CnoidalWaveVelY_ = 8
integer, parameter :: SolitartyWaveVelY_ = 9
!wave cell types
integer, parameter :: EnteringWaveCell_ = 1
integer, parameter :: LeavingWaveCell_ = 2
integer, parameter :: InteriorWaveCell_ = 3
integer, parameter :: ExteriorWaveCell_ = 4
!Parameter-----------------------------------------------------------------
character(LEN = StringLength), parameter :: BeginProfile = '<BeginProfile>'
character(LEN = StringLength), parameter :: EndProfile = '<EndProfile>'
character(LEN = StringLength), parameter :: BeginLayers = '<BeginLayers>'
character(LEN = StringLength), parameter :: EndLayers = '<EndLayers>'
character(LEN = StringLength), parameter :: BeginDepth = '<BeginDepth>'
character(LEN = StringLength), parameter :: EndDepth = '<EndDepth>'
character(LEN = StringLength), parameter :: BeginTimes = '<BeginTime>'
character(LEN = StringLength), parameter :: EndTimes = '<EndTime>'
character(LEN = StringLength), parameter :: BeginProfileValues= '<BeginProfileValues>'
character(LEN = StringLength), parameter :: EndProfileValues = '<EndProfileValues>'
character(LEN = StringLength), parameter :: BeginMTSBlock = '<BeginStation>'
character(LEN = StringLength), parameter :: EndMTSBlock = '<EndStation>'
character(LEN = StringLength), parameter :: Begin_files = '<<BeginFiles>>'
character(LEN = StringLength), parameter :: End_files = '<<EndFiles>>'
character(LEN = StringLength), parameter :: Begin_files_2 = '<<<BeginFiles>>>'
character(LEN = StringLength), parameter :: End_files_2 = '<<<EndFiles>>>'
!Types---------------------------------------------------------------------
type T_Layers
real, dimension(:), pointer :: Values => null()
end type T_Layers
type T_Boxes
character(PathLength) :: FileName = null_str
integer :: ObjBoxDif = null_int
real, dimension(:), pointer :: Values => null()
end type T_Boxes
type T_ASCIIFile
character(PathLength) :: FileName = null_str
integer :: GridDataID = null_int
type (T_ASCIIFile), pointer :: Next => null()
type (T_ASCIIFile), pointer :: Prev => null()
end type T_ASCIIFile
type T_Sponge
real :: OutValue = null_real
integer :: Cells = null_int
logical :: Growing = .false.
integer :: Evolution = null_int
!1 - South; 2 - North; 3 - West; 4 - East
logical, dimension(1:4) :: OpenBordersON = .true.
end type T_Sponge
type T_EnteringCell
integer, dimension(:), pointer :: i => null()
integer, dimension(:), pointer :: j => null()
real, dimension(:,:), pointer :: dx => null()
real, dimension(:,:), pointer :: dy => null()
real, dimension(:), pointer :: TimeLag => null()
integer :: nCells = null_int
integer :: iStart = null_int
integer :: jStart = null_int
integer :: n = null_int
end type T_EnteringCell
type T_AnalyticWave
logical :: ON = .false.
real :: Amplitude = null_real
real :: Direction = null_real
real :: Period = null_real
real :: AverageValue = null_real
real :: DepthValue = null_real
real :: CoefValue = null_real
logical :: SlowStartON = .false.
! WaveType = 1 (Sine), WaveType = 2 (Cnoidal), WaveType = 3 (solitary)
integer :: WaveType = null_int
real, dimension(:,:), pointer :: X2D => null()
real, dimension(:,:), pointer :: Celerity => null()
real, dimension(:,:), pointer :: AmpAux => null()
! EnteringWaveCell_ = 1
! LeavingWaveCell_ = 2
! InteriorWaveCell_ = 3
! ExteriorWaveCell_ = 4
integer, dimension(:,:), pointer :: CellType => null()
integer, dimension(:,:), pointer :: TlagMissing => null()
type (T_EnteringCell) :: EnteringCell
end type T_AnalyticWave
type T_TimeSerie
character(PathLength) :: FileName = null_str
integer :: ObjTimeSerie = 0
integer :: Column = null_int
type (T_Time) :: NextTime, &
PreviousTime
integer :: NextInstant = 0, &
PreviousInstant = 0
real :: PreviousValue = 0., &
NextValue = 0., &
CurrentValue = 0.
integer :: NumberOfInstants = 0
logical :: RemainsConstant = .false.
real :: PredictedDT = -null_real
real :: DTForNextEvent = -null_real
real :: DTForNextDataset = -null_real
type(T_Time) :: NextEventStart
type(T_Time) :: NextEventEnd
type(T_Time) :: TimeOfNextDataset
integer :: InstantOfNextDataset = null_int
real :: NextValueForDTPred = -null_real
type (T_TimeSerie), pointer :: Next => null()
type (T_TimeSerie), pointer :: Prev => null()
end type T_TimeSerie
type T_ProfileTimeSerie
character(PathLength) :: FileName = null_str
type (T_Time) :: NextTime, PreviousTime
integer :: NextInstant = null_int, &
PreviousInstant = null_int
real, dimension(:,:,:), pointer :: PreviousField3D => null(), &
NextField3D => null()
real, dimension(:,: ), pointer :: Values => null(), &
Depths => null()
type(T_Time), dimension(: ), pointer :: TimeInstants => null()
integer :: NumberOfInstants = null_int, &
nValues = null_int, &
nDepths = null_int, &
FirstInstant = null_int, &
LastInstant = null_int
logical :: CyclicTimeON = .false.
end type T_ProfileTimeSerie
type T_Station
character(PathLength) :: FileName = null_str
integer :: ObjTimeSerie = 0
integer :: Column = null_int
integer :: FillID = null_int
logical :: RemainConstant = .false.
logical :: ValueIsDefined = .false.
real :: NewValue = -null_real
integer :: NumberOfInstants = 0
type(T_Time) :: NextEventStart
type(T_Time) :: NextEventEnd
type(T_Time) :: NextTime
type(T_Time) :: PreviousTime
integer :: NextInstant = 0
integer :: PreviousInstant = 0
real :: PredictedDT = -null_real
real :: DTForNextEvent = 0.
real :: DTForNextDataset = -null_real
real :: NextValue = 0.
real :: PreviousValue = 0.
real :: NextValueForDTPred = 0.
end type T_Station
type T_MultiTimeSerie
integer, dimension(:,:), allocatable :: FillGrid2D
integer, dimension(:,:,:), allocatable :: FillGrid3D
integer :: DataProcessing = null_int
type(T_Station), dimension(:), allocatable :: StationsList
integer :: NumberOfSources = 0
end type T_MultiTimeserie
!Generic 4D
type T_Generic4D
logical :: ON = .false.
logical :: ReadFromTimeSerie = .false.
integer :: ObjTimeSerie = null_int
integer :: TimeSerieColumn = null_int
real :: CurrentValue = null_real
end type
type T_Field4D
character(PathLength) :: FileName = null_str, &
VGroupPath = null_str, &
FieldName = null_str
real :: MultiplyingFactor = null_real
logical :: HasMultiplyingFactor = .false.
real :: AddingFactor = null_real
logical :: HasAddingFactor = .false.
type (T_Time) :: NextTime, PreviousTime
type (T_Time) :: StartTime, EndTime
real :: Next4DValue = FillValueReal
real :: Previous4DValue = FillValueReal
integer :: NextInstant = null_int, &
PreviousInstant = null_int
real, dimension(:,: ), pointer :: PreviousField2D => null(), &
NextField2D => null(), &
Array2D => null()
real, dimension(:,:,:), pointer :: PreviousField3D => null(), &
NextField3D => null(), &
ReadField3D => null(), &
Array3D => null()
real :: PredictedDT = -null_real
real :: DTForNextEvent = -null_real
real :: DTForNextDataset = -null_real
type(T_Time) :: NextEventStart
type(T_Time) :: NextEventEnd
type(T_Time) :: TimeOfNextDataset
integer :: InstantOfNextDataset = null_int
real :: NextValueForDTPred = -null_real
integer :: ObjHDF5 = 0
logical :: RemainsConstant = .false.
integer :: NumberOfInstants = null_int
logical :: CyclicTimeON = .false.
logical :: From2Dto3D = .false.
logical :: From3Dto2D = .false.
type(T_Generic4D) :: Generic4D
!logical :: ArgumentFileName = .false.
integer :: ObjField4D = 0
logical :: Field4D = .false.
logical :: HarmonicsON = .false.
real :: HarmonicsDT = null_real
logical :: SpatialInterpolON = .false.
logical :: InterpolOnlyVertically = .false.
logical :: GenericYear = .false.
integer :: Ncells
real, dimension(:), pointer :: X => null()
real, dimension(:), pointer :: Y => null()
real, dimension(:), pointer :: Z => null()
real, dimension(:), pointer :: Prop => null()
logical, dimension(:), pointer :: NoData => null()
logical :: Extrapolate = .false.
integer :: ExtrapolateMethod = null_int
character(len=PathLength), dimension(:), pointer :: FileNameList => null()
type (T_Field4D), pointer :: Next => null()
type (T_Field4D), pointer :: Prev => null()
end type T_Field4D
private :: T_FillMatrix
type T_FillMatrix
integer :: InstanceID
type (T_Size2D) :: Size2D, WorkSize2D
type (T_Size3D) :: Size3D, WorkSize3D
type (T_PropertyID) :: PropertyID
integer :: Dim = null_int
integer :: TypeZUV = null_int
integer :: TimeEvolution = null_int
integer :: SpaceEvolution = null_int
integer :: InitializationMethod = null_int
integer :: InitializationDefault = null_int
logical :: RemainsConstant = .false.
logical :: VectorialProp = .false. ! w/ x and y fields (orig & rotated)
logical :: RotateAngleToGrid = .false. !scalar w/ original and rotated field
! logical :: AccumulatedValue = .false.
! logical :: NoInterpol = .false.
! integer :: ValuesType
logical :: InterpolateValues = .false.
logical :: AccumulateValues = .false.
logical :: UseOriginalValues = .false.
logical :: PreviousInstantValues = .false.
logical :: IgnoreNoDataPoint = .false.
integer :: PredictDTMethod !1 for old method, 2 for new method (for rain, mainly)
real :: NoDataValue = null_real
logical :: Backtracking = .false.
character(len=StringLength) :: OverrideValueKeyword = null_str
logical :: OverrideValueKeywordON = .false.
real :: MinForDTDecrease = AllmostZero
!real :: DefaultValue = null_real
real, dimension(3) :: DefaultValue = null_real
real :: PredictedDT = -null_real
real :: DTForNextEvent = -null_real
real :: DTForNextDataset = -null_real
!type(T_Time) :: NextEventStart
!type(T_Time) :: NextEventEnd
!type(T_Time) :: TimeOfNextDataset
!integer :: InstantOfNextDataset
logical :: ValueIsUsedForDTPrediction = .false.
real :: NextValueForDTPred
real, dimension(:, : ), pointer :: Matrix2D => null()
real, dimension(:, : ), pointer :: Matrix2DX => null() !input vectorial field x (meridian value)
real, dimension(:, : ), pointer :: Matrix2DY => null() !input vectorial field y (parallel value)
real, dimension(:, : ), pointer :: Matrix2DU => null() !model grid vectorial field U (cell ref)
real, dimension(:, : ), pointer :: Matrix2DV => null() !model grid vectorial field V (cell ref)
real, dimension(:, : ), pointer :: Matrix2DFieldAngle => null() !input angle field in nautic or current ref
real, dimension(:, : ), pointer :: Matrix2DCellAngle => null() !cell angle field in cell ref
real, dimension(:, :, :), pointer :: Matrix3D => null()
real, dimension(:, :, :), pointer :: Matrix3DX => null()
real, dimension(:, :, :), pointer :: Matrix3DY => null()
real, dimension(:, :, :), pointer :: Matrix3DZ => null()
real, dimension(:, :, :), pointer :: Matrix3DU => null()
real, dimension(:, :, :), pointer :: Matrix3DV => null()
real, dimension(:, :, :), pointer :: Matrix3DW => null()
real, dimension(:, :, :), pointer :: Matrix3DFieldAngle => null()
real, dimension(:, :, :), pointer :: Matrix3DCellAngle => null()
integer, dimension(:, : ), pointer :: PointsToFill2D => null()
integer, dimension(:, :, :), pointer :: PointsToFill3D => null()
logical :: UseZ = .false. !In case of 3D,use Z component?
type(T_Time) :: BeginTime, EndTime
logical :: ArgumentFileName = .false.
character(len=StringLength), dimension(2) :: FileNameHDF = null_str
type(T_ASCIIFile), pointer :: FirstASCIIFile
type(T_TimeSerie), pointer :: FirstTimeSerie
type(T_Field4D), pointer :: FirstHDF
integer :: nASCIIFiles = null_int
integer :: nTimeSeries = null_int
integer :: nHDFs = null_int
logical :: CheckDates = .true.
character(len=PathLength) :: SpongeFILE_DT = null_str
!Initialization Methods
type (T_Layers ) :: Layers
type (T_Boxes ) :: Boxes
type (T_TimeSerie) :: TimeSerie
type (T_ASCIIFile) :: ASCIIFile
type (T_Sponge ) :: Sponge
type (T_AnalyticWave) :: AnalyticWave
type (T_Field4D ) :: HDF
type (T_ProfileTimeSerie) :: ProfileTimeSerie
type (T_MultiTimeSerie) :: MultiTimeSerie
integer :: ObjEnterData = 0
integer :: ObjTime = 0
integer :: ObjHorizontalGrid = 0
integer :: ObjGeometry = 0
type(T_FillMatrix), pointer :: Next => null()
end type T_FillMatrix
!Global Module Variables
type (T_FillMatrix), pointer :: FirstObjFillMatrix => null()
type (T_FillMatrix), pointer :: Me => null()
!--------------------------------------------------------------------------
contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONS
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine ConstructFillMatrix2D(PropertyID, EnterDataID, TimeID, &
HorizontalGridID, ExtractType, PointsToFill2D, &
Matrix2D, TypeZUV, Matrix2DInputRef, FileNameHDF, &
ObjFillMatrix, OverrideValueKeyword, ClientID, &
PredictDTMethod, MinForDTDecrease, &
ValueIsUsedForDTPrediction, CheckDates, &
RotateAngleToGrid, SpongeFILE_DT, STAT)
!Arguments---------------------------------------------------------------
integer :: EnterDataID
integer :: TimeID
integer :: HorizontalGridID
integer :: ExtractType
integer, dimension(:, :), pointer :: PointsToFill2D
real, dimension(:, :), pointer :: Matrix2D
real, dimension(:, :), pointer, optional :: Matrix2DInputRef !original field (e.g. angle)
integer :: TypeZUV
type (T_PropertyID) :: PropertyID
integer, optional, intent(IN) :: ClientID
character(*), optional, intent(IN ) :: FileNameHDF, OverrideValueKeyword
integer, optional, intent(INOUT) :: ObjFillMatrix
logical, optional, intent(IN) :: CheckDates
integer, optional, intent(OUT) :: STAT
integer, optional, intent(IN ) :: PredictDTMethod
real, optional, intent(IN ) :: MinForDTDecrease
logical, optional, intent(IN ) :: ValueIsUsedForDTPrediction
logical, optional, intent(IN ) :: RotateAngleToGrid
character(*), optional, intent(IN ) :: SpongeFILE_DT
!Local-------------------------------------------------------------------
integer :: ready_, STAT_, STAT_CALL, nUsers, ObjFillMatrix_
integer :: PredictDTMethod_, Referential
type (T_PropertyID), pointer :: Prop
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mFillMatrix_)) then
nullify (FirstObjFillMatrix)
call RegisterModule (mFillMatrix_)
endif
if (present(ObjFillMatrix)) then
ObjFillMatrix_ = ObjFillMatrix
else
ObjFillMatrix_ = PropertyID%ObjFillMatrix
endif
call Ready(ObjFillMatrix_, ready_)
cd0 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
if (present(CheckDates)) then
Me%CheckDates = CheckDates
endif
if (present(SpongeFILE_DT)) then
Me%SpongeFILE_DT = trim(SpongeFILE_DT)
else
Me%SpongeFILE_DT = null_str
endif
!~ if (Check_Vectorial_Property(PropertyID%IDNumber)) then
if (PropertyID%IsVectorial) then
write(*,*) 'Constructing vectorial property but expected scalar'
stop 'ConstructFillMatrix2D - ModuleFillMatrix - ERR00'
endif
if (present(PredictDTMethod)) then
PredictDTMethod_ = PredictDTMethod
else
PredictDTMethod_ = 1
endif
if (present(MinForDTDecrease)) then
Me%MinForDTDecrease = MinForDTDecrease
endif
if (present(ValueIsUsedForDTPrediction)) then
Me%ValueIsUsedForDTPrediction = ValueIsUsedForDTPrediction
else
Me%ValueIsUsedForDTPrediction = .false.
endif
Me%ObjEnterData = AssociateInstance (mENTERDATA_, EnterDataID )
Me%ObjTime = AssociateInstance (mTIME_, TimeID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)
! JPW 2012-01-28: Use HorizontalGridSize to set size of the matrix.
! Using ubound may cause inconsistency with padded matrices (making IUB larger than the actual grid bound).
! PointsToFill2D should have the same dimensions as HorizontalGrid anyway.
! Matrices that are used for the Thomas algorithm are padded if _USE_PAGELOCKED (CUDA only) or _PAD_MATRICES (Fortran) is defined
call GetHorizontalGridSize(HorizontalGridID, Me%Size2D, Me%WorkSize2D, STAT = STAT_)
! Me%Size2D%ILB = lbound(PointsToFill2D, dim = 1)
! Me%Size2D%IUB = ubound(PointsToFill2D, dim = 1)
! Me%Size2D%JLB = lbound(PointsToFill2D, dim = 2)
! Me%Size2D%JUB = ubound(PointsToFill2D, dim = 2)
! Me%WorkSize2D%ILB = Me%Size2D%ILB + 1
! Me%WorkSize2D%IUB = Me%Size2D%IUB - 1
! Me%WorkSize2D%JLB = Me%Size2D%JLB + 1
! Me%WorkSize2D%JUB = Me%Size2D%JUB - 1
Me%Size3D = T_Size3D(null_int, null_int, null_int, null_int, null_int, null_int)
Me%Dim = Dim2D
Me%TypeZUV = TypeZUV
Me%PropertyID = PropertyID
if (present(RotateAngleToGrid)) then
Me%RotateAngleToGrid = RotateAngleToGrid
else
if (Me%PropertyID%IsAngle) then
Me%RotateAngleToGrid = .true.
endif
endif
Me%Matrix2D => Matrix2D
Me%PointsToFill2D => PointsToFill2D
!!get the orginal field. will be given to user to output
!~ if (Check_Angle_Property(Me%PropertyID%IDNumber)) then
if (Me%RotateAngleToGrid) then
if (.not. present(Matrix2DInputRef)) then
write(*,*) 'Constructing angle property but not given original field'
stop 'ConstructFillMatrix2D - ModuleFillMatrix - ERR10'
endif
!point to input angle and cell angle (rotated to cell ref)
Me%Matrix2DFieldAngle => Matrix2DInputRef
Me%Matrix2DCellAngle => Matrix2D
where (PointsToFill2D == WaterPoint) Me%Matrix2DFieldAngle = null_real
where (PointsToFill2D == WaterPoint) Me%Matrix2DCellAngle = null_real
endif
where (PointsToFill2D == WaterPoint) Me%Matrix2D = null_real
if (Me%TypeZUV == TypeU_) then
Me%Size2D%JUB = Me%Size2D%JUB + 1
Me%WorkSize2D%JUB = Me%WorkSize2D%JUB + 1
endif
if (Me%TypeZUV == TypeV_) then
Me%Size2D%IUB = Me%Size2D%IUB + 1
Me%WorkSize2D%IUB = Me%WorkSize2D%IUB + 1
endif
if (present(FileNameHDF)) then
Me%ArgumentFileName = .true.
Me%FileNameHDF(1) = trim(FileNameHDF)
else
Me%ArgumentFileName = .false.
endif
if(present(OverrideValueKeyword))then
Me%OverrideValueKeyword = trim(adjustl(OverrideValueKeyword))
Me%OverrideValueKeywordON = .true.
else
Me%OverrideValueKeywordON = .false.
end if
if (present(ClientID)) then
call ReadOptions (ExtractType, &
PointsToFill2D = PointsToFill2D, &
PredictDTMethod = PredictDTMethod_, &
ClientID = ClientID)
else
call ReadOptions (ExtractType, &
PointsToFill2D = PointsToFill2D, &
PredictDTMethod = PredictDTMethod_)
endif
!Is this a angle property? convert to cell referential angle
if (Me%RotateAngleToGrid) then
!angle referential
Prop => Me%PropertyID
Referential = Get_Angle_Referential(Prop)
!!Need to rotate input field
call RotateAngleFieldToGrid(HorizontalGridID = Me%ObjHorizontalGrid, &
AngleIn = Me%Matrix2DFieldAngle, &
InReferential = Referential, &
AngleOut = Me%Matrix2DCellAngle, &
WaterPoints2D = PointsToFill2D, &
Rotate = .true., &
STAT = STAT_CALL)
endif
nUsers = DeassociateInstance(mENTERDATA_, Me%ObjEnterData)
if (nUsers == 0) stop 'ConstructFillMatrix2D - ModuleFillMatrix - ERR20'
if(Me%TimeEvolution == None)then
PropertyID%SolutionFromFile = .false.
else
PropertyID%SolutionFromFile = .true.
end if
!Returns ID
ObjFillMatrix_ = Me%InstanceID
if (present(ObjFillMatrix)) then
ObjFillMatrix = ObjFillMatrix_
else
PropertyID%ObjFillMatrix = ObjFillMatrix_
endif
nullify(Me%Matrix2D )
nullify(Me%Matrix2DFieldAngle )
nullify(Me%Matrix2DCellAngle )
nullify(Me%PointsToFill2D)
STAT_ = SUCCESS_
else cd0
stop 'ConstructFillMatrix2D - ModuleFillMatrix - ERR02'
end if cd0
if (present(STAT)) STAT = STAT_
!----------------------------------------------------------------------
end subroutine ConstructFillMatrix2D
!----------------------------------------------------------------------
subroutine ConstructFillMatrix2DVectorial(PropertyID, EnterDataID, TimeID, &
HorizontalGridID, ExtractType, PointsToFill2D, &
Matrix2DU, Matrix2DV, Matrix2DX, Matrix2DY, &
TypeZUV, FileNameHDF, ObjFillMatrix, &
OverrideValueKeyword, ClientID, PredictDTMethod, &
MinForDTDecrease, ValueIsUsedForDTPrediction, CHeckDates, STAT)
!Arguments---------------------------------------------------------------
integer :: EnterDataID
integer :: TimeID
integer :: HorizontalGridID
integer :: ExtractType
integer, dimension(:, :), pointer :: PointsToFill2D
real, dimension(:, :), pointer :: Matrix2DU
real, dimension(:, :), pointer :: Matrix2DV
real, dimension(:, :), pointer :: Matrix2DX
real, dimension(:, :), pointer :: Matrix2DY
integer :: TypeZUV
type (T_PropertyID) :: PropertyID
integer, optional, intent(IN) :: ClientID
character(*), dimension(2), optional, intent(IN ) :: FileNameHDF
character(*), optional, intent(IN ) :: OverrideValueKeyword
integer, optional, intent(INOUT) :: ObjFillMatrix
logical, optional, intent(IN) :: CheckDates
integer, optional, intent(OUT) :: STAT
integer, optional, intent(IN ) :: PredictDTMethod
real, optional, intent(IN ) :: MinForDTDecrease
logical, optional, intent(IN ) :: ValueIsUsedForDTPrediction
!Local-------------------------------------------------------------------
integer :: ready_, STAT_, STAT_CALL, nUsers, ObjFillMatrix_
integer :: PredictDTMethod_
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mFillMatrix_)) then
nullify (FirstObjFillMatrix)
call RegisterModule (mFillMatrix_)
endif
if (present(ObjFillMatrix)) then
ObjFillMatrix_ = ObjFillMatrix
else
ObjFillMatrix_ = PropertyID%ObjFillMatrix
endif
call Ready(ObjFillMatrix_, ready_)
cd0 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
if (present(CheckDates)) then
Me%CheckDates = CheckDates
endif
Me%VectorialProp = .true.
!~ if (.not. Check_Vectorial_Property(PropertyID%IDNumber)) then
if (.not. PropertyID%IsVectorial) then
write(*,*) 'Constructing scalar property but expected vectorial'
stop 'ConstructFillMatrix2DVectorial - ModuleFillMatrix - ERR00'
endif
if (present(PredictDTMethod)) then
PredictDTMethod_ = PredictDTMethod
else
PredictDTMethod_ = 1
endif
if (present(MinForDTDecrease)) then
Me%MinForDTDecrease = MinForDTDecrease
endif
if (present(ValueIsUsedForDTPrediction)) then
Me%ValueIsUsedForDTPrediction = ValueIsUsedForDTPrediction
else
Me%ValueIsUsedForDTPrediction = .false.
endif
Me%ObjEnterData = AssociateInstance (mENTERDATA_, EnterDataID )
Me%ObjTime = AssociateInstance (mTIME_, TimeID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)
! JPW 2012-01-28: Use HorizontalGridSize to set size of the matrix.
! Using ubound may cause inconsistency with padded matrices (making IUB larger than the actual grid bound).
! PointsToFill2D should have the same dimensions as HorizontalGrid anyway.
! Matrices that are used for the Thomas algorithm are padded if _USE_PAGELOCKED (CUDA only) or _PAD_MATRICES (Fortran) is defined
call GetHorizontalGridSize(HorizontalGridID, Me%Size2D, Me%WorkSize2D, STAT = STAT_)
! Me%Size2D%ILB = lbound(PointsToFill2D, dim = 1)
! Me%Size2D%IUB = ubound(PointsToFill2D, dim = 1)
! Me%Size2D%JLB = lbound(PointsToFill2D, dim = 2)
! Me%Size2D%JUB = ubound(PointsToFill2D, dim = 2)
! Me%WorkSize2D%ILB = Me%Size2D%ILB + 1
! Me%WorkSize2D%IUB = Me%Size2D%IUB - 1
! Me%WorkSize2D%JLB = Me%Size2D%JLB + 1
! Me%WorkSize2D%JUB = Me%Size2D%JUB - 1
Me%Size3D = T_Size3D(null_int, null_int, null_int, null_int, null_int, null_int)
Me%Dim = Dim2D
Me%TypeZUV = TypeZUV
Me%PropertyID = PropertyID
!Point result matrixes to argument (user wants rotated to grid result)
Me%Matrix2DU => Matrix2DU
Me%Matrix2DV => Matrix2DV
Me%PointsToFill2D => PointsToFill2D
!Point input matrixes to argument (user inpt field)
Me%Matrix2DX => Matrix2DX
Me%Matrix2DY => Matrix2DY
!Start matrice
!where (PointsToFill2D == WaterPoint) Me%Matrix2D = null_real
where (PointsToFill2D == WaterPoint) Me%Matrix2DU = null_real
where (PointsToFill2D == WaterPoint) Me%Matrix2DV = null_real
where (PointsToFill2D == WaterPoint) Me%Matrix2DX = null_real
where (PointsToFill2D == WaterPoint) Me%Matrix2DY = null_real
if (Me%TypeZUV == TypeU_) then
Me%Size2D%JUB = Me%Size2D%JUB + 1
Me%WorkSize2D%JUB = Me%WorkSize2D%JUB + 1
endif
if (Me%TypeZUV == TypeV_) then
Me%Size2D%IUB = Me%Size2D%IUB + 1
Me%WorkSize2D%IUB = Me%WorkSize2D%IUB + 1
endif
if (present(FileNameHDF)) then
Me%ArgumentFileName = .true.
Me%FileNameHDF(1) = trim(FileNameHDF(1))
Me%FileNameHDF(2) = trim(FileNameHDF(2))
else
Me%ArgumentFileName = .false.
endif
if(present(OverrideValueKeyword))then
Me%OverrideValueKeyword = trim(adjustl(OverrideValueKeyword))
Me%OverrideValueKeywordON = .true.
else
Me%OverrideValueKeywordON = .false.
end if
!Here the Me%Matrix2DX and Me%Matrix2DY from user input field are filled
if (present(ClientID)) then
call ReadOptions (ExtractType, &
PointsToFill2D = PointsToFill2D, &
PredictDTMethod = PredictDTMethod_, &
ClientID = ClientID)
else
call ReadOptions (ExtractType, &
PointsToFill2D = PointsToFill2D, &
PredictDTMethod = PredictDTMethod_)
endif
nUsers = DeassociateInstance(mENTERDATA_, Me%ObjEnterData)
if (nUsers == 0) stop 'ConstructFillMatrix2D - ModuleFillMatrix - ERR20'
!!Need to rotate input field (Me%Matrix2DX and Me%Matrix2DY) to grid (Me%Matrix2DU and Me%Matrix2DV)
call RotateVectorFieldToGrid(HorizontalGridID = Me%ObjHorizontalGrid, &
VectorInX = Me%Matrix2DX, &
VectorInY = Me%Matrix2DY, &
VectorOutX = Me%Matrix2DU, &
VectorOutY = Me%Matrix2DV, &
WaterPoints2D = PointsToFill2D, &
RotateX = .true., &
RotateY = .true., &
STAT = STAT_CALL)
if(Me%TimeEvolution == None)then
PropertyID%SolutionFromFile = .false.
else
PropertyID%SolutionFromFile = .true.
end if
!Returns ID
ObjFillMatrix_ = Me%InstanceID
if (present(ObjFillMatrix)) then
ObjFillMatrix = ObjFillMatrix_
else
PropertyID%ObjFillMatrix = ObjFillMatrix_
endif
!nullify will remove pointer to Matrix2DU, Matrix2DV and X and Y (Module fields) since they were already modifed
nullify(Me%Matrix2D )
nullify(Me%Matrix2DU )
nullify(Me%Matrix2DV )
nullify(Me%Matrix2DX )
nullify(Me%Matrix2DY )
nullify(Me%PointsToFill2D)
STAT_ = SUCCESS_
else cd0
stop 'ConstructFillMatrix2D - ModuleFillMatrix - ERR02'
end if cd0
if (present(STAT)) STAT = STAT_
!----------------------------------------------------------------------
end subroutine ConstructFillMatrix2DVectorial
!----------------------------------------------------------------------
subroutine ConstructFillMatrix3D(PropertyID, EnterDataID, TimeID, HorizontalGridID, &
GeometryID, ExtractType, PointsToFill3D, Matrix3D, &
TypeZUV, Matrix3DInputRef, FillMatrix, FileNameHDF,&
ObjFillMatrix, OverrideValueKeyword, ClientID, &
predictDTMethod, MinForDTDecrease, &
ValueIsUsedForDTPrediction, CheckDates, &
RotateAngleToGrid, SpongeFILE_DT, STAT)
!Arguments---------------------------------------------------------------
type (T_PropertyID) :: PropertyID
integer :: EnterDataID
integer :: TimeID
integer :: HorizontalGridID
integer :: GeometryID
integer :: ExtractType
integer, dimension(:, :, :), pointer :: PointsToFill3D
real, dimension(:, :, :), pointer :: Matrix3D
real, dimension(:, :, :), pointer, optional :: Matrix3DInputRef !original field (e.g. angle)
integer :: TypeZUV
integer, optional, intent(IN) :: ClientID
real , optional, intent(IN ) :: FillMatrix
character(*), optional, intent(IN ) :: FileNameHDF, OverrideValueKeyword
integer, optional, intent(INOUT) :: ObjFillMatrix
logical, optional, intent(IN) :: CheckDates
integer, optional, intent(OUT) :: STAT
integer, optional, intent(IN ) :: PredictDTMethod
real, optional, intent(IN ) :: MinForDTDecrease
logical, optional, intent(IN ) :: ValueIsUsedForDTPrediction
logical, optional, intent(IN ) :: RotateAngleToGrid
character(*), optional, intent(IN ) :: SpongeFILE_DT
!Local-------------------------------------------------------------------
real :: FillMatrix_
integer :: ready_, STAT_, STAT_CALL, nUsers, ObjFillMatrix_
integer :: PredictDTMethod_, Referential
type (T_PropertyID), pointer :: Prop
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mFillMatrix_)) then
nullify (FirstObjFillMatrix)
call RegisterModule (mFillMatrix_)
endif
if (present(ObjFillMatrix)) then
ObjFillMatrix_ = ObjFillMatrix
else
ObjFillMatrix_ = PropertyID%ObjFillMatrix
endif
call Ready(ObjFillMatrix_, ready_)
cd0 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
if (present(RotateAngleToGrid)) then
Me%RotateAngleToGrid = RotateAngleToGrid
else
if (Me%PropertyID%IsAngle) then
Me%RotateAngleToGrid = .true.
endif
endif
if (present(CheckDates)) then
Me%CheckDates = CheckDates
endif
if (present(SpongeFILE_DT)) then
Me%SpongeFILE_DT = trim(SpongeFILE_DT)
else
Me%SpongeFILE_DT = null_str
endif
!~ if (Check_Vectorial_Property(PropertyID%IDNumber)) then
if (PropertyID%IsVectorial) then
write(*,*) 'Constructing vectorial property but expected scalar'
stop 'ConstructFillMatrix3D - ModuleFillMatrix - ERR00'
endif
if (present(PredictDTMethod)) then
PredictDTMethod_ = PredictDTMethod
else
PredictDTMethod_ = 1
endif
if (present(MinForDTDecrease)) then
Me%MinForDTDecrease = MinForDTDecrease
endif
if (present(ValueIsUsedForDTPrediction)) then
Me%ValueIsUsedForDTPrediction = ValueIsUsedForDTPrediction
else
Me%ValueIsUsedForDTPrediction = .false.
endif
Me%ObjEnterData = AssociateInstance (mENTERDATA_, EnterDataID )
Me%ObjTime = AssociateInstance (mTIME_, TimeID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)
Me%ObjGeometry = AssociateInstance (mGEOMETRY_, GeometryID )
! JPW 2012-01-28: Use GeometrySize to set size of the matrix.
! Using ubound may cause inconsistency with padded matrices (making IUB larger than the actual grid bound).
! PointsToFill3D should have the same dimensions as Geometry anyway.
! Matrices that are used for the Thomas algorithm are padded if _USE_PAGELOCKED (CUDA only) or _PAD_MATRICES (Fortran) is defined
call GetGeometrySize(GeometryID , Me%Size3D, Me%WorkSize3D, STAT_)
! Me%Size3D%ILB = lbound(PointsToFill3D, dim = 1)
! Me%Size3D%IUB = ubound(PointsToFill3D, dim = 1)
! Me%Size3D%JLB = lbound(PointsToFill3D, dim = 2)
! Me%Size3D%JUB = ubound(PointsToFill3D, dim = 2)
! Me%Size3D%KLB = lbound(PointsToFill3D, dim = 3)
! Me%Size3D%KUB = ubound(PointsToFill3D, dim = 3)
! Me%WorkSize3D%ILB = Me%Size3D%ILB + 1
! Me%WorkSize3D%IUB = Me%Size3D%IUB - 1
! Me%WorkSize3D%JLB = Me%Size3D%JLB + 1
! Me%WorkSize3D%JUB = Me%Size3D%JUB - 1
! Me%WorkSize3D%KLB = Me%Size3D%KLB + 1
! Me%WorkSize3D%KUB = Me%Size3D%KUB - 1
Me%Size2D = T_Size2D(null_int, null_int, null_int, null_int)
Me%Dim = Dim3D
Me%TypeZUV = TypeZUV
Me%PropertyID = PropertyID
Me%Matrix3D => Matrix3D
Me%PointsToFill3D => PointsToFill3D
!!get the orginal field. will be given to user to output
!~ if (Check_Angle_Property(Me%PropertyID%IDNumber)) then
if (Me%RotateAngleToGrid) then
if (.not. present(Matrix3DInputRef)) then
write(*,*) 'Constructing angle property but not given original field'
stop 'ConstructFillMatrix3D - ModuleFillMatrix - ERR10'
endif
!point to input angle and cell angle (rotated to cell ref)
Me%Matrix3DFieldAngle => Matrix3DInputRef
Me%Matrix3DCellAngle => Matrix3D
where (PointsToFill3D == WaterPoint) Me%Matrix3DFieldAngle = null_real
where (PointsToFill3D == WaterPoint) Me%Matrix3DCellAngle = null_real
endif
if (present(FillMatrix)) then
FillMatrix_ = FillMatrix
else
FillMatrix_ = null_real
endif
where (PointsToFill3D == WaterPoint) Me%Matrix3D = FillMatrix_
if (Me%TypeZUV == TypeU_) then
Me%Size3D%JUB = Me%Size3D%JUB + 1
Me%WorkSize3D%JUB = Me%WorkSize3D%JUB + 1
endif
if (Me%TypeZUV == TypeV_) then
Me%Size3D%IUB = Me%Size3D%IUB + 1
Me%WorkSize3D%IUB = Me%WorkSize3D%IUB + 1
endif
!Specific of the Vertical_Z matrix see ModuleGeometry
if (Me%TypeZUV == TypeW_) then
Me%WorkSize3D%KLB = Me%WorkSize3D%KLB - 1
endif
if (present(FileNameHDF)) then
Me%ArgumentFileName = .true.
Me%FileNameHDF(1) = trim(FileNameHDF)
else
Me%ArgumentFileName = .false.
endif
if(present(OverrideValueKeyword))then
Me%OverrideValueKeyword = trim(adjustl(OverrideValueKeyword))
Me%OverrideValueKeywordON = .true.
else
Me%OverrideValueKeywordON = .false.
end if
if (present(ClientID)) then
call ReadOptions (ExtractType, &
PointsToFill3D = PointsToFill3D, &
PredictDTMethod = PredictDTMethod_, &
ClientID = ClientID)
else
call ReadOptions (ExtractType, &
PointsToFill3D = PointsToFill3D, &
PredictDTMethod = PredictDTMethod_)
endif
!Is this a angle property? convert to cell referential angle
if (Me%RotateAngleToGrid) then
!angle referential
Prop => Me%PropertyID
Referential = Get_Angle_Referential(Prop)
!!Need to rotate input field
call RotateAngleFieldToGrid(HorizontalGridID = Me%ObjHorizontalGrid, &
AngleIn = Me%Matrix3DFieldAngle, &
InReferential = Referential, &
AngleOut = Me%Matrix3DCellAngle, &
WaterPoints3D = PointsToFill3D, &
Rotate = .true., &
KLB = Me%WorkSize3D%KLB, &
KUB = Me%WorkSize3D%KUB, &
STAT = STAT_CALL)
endif
nUsers = DeassociateInstance(mENTERDATA_, Me%ObjEnterData)
if (nUsers == 0) stop 'ConstructFillMatrix3D - ModuleFillMatrix - ERR01'
if(Me%TimeEvolution == None)then
PropertyID%SolutionFromFile = .false.
else
PropertyID%SolutionFromFile = .true.
end if
!Returns ID
ObjFillMatrix_ = Me%InstanceID
if (present(ObjFillMatrix)) then
ObjFillMatrix = ObjFillMatrix_
else
PropertyID%ObjFillMatrix = ObjFillMatrix_
endif
nullify(Me%Matrix3D )
nullify(Me%Matrix3DFieldAngle )
nullify(Me%Matrix3DCellAngle )
nullify(Me%PointsToFill3D)
STAT_ = SUCCESS_
else cd0
stop 'ConstructFillMatrix3D - ModuleFillMatrix - ERR02'
end if cd0
if (present(STAT)) STAT = STAT_
!----------------------------------------------------------------------
end subroutine ConstructFillMatrix3D
!--------------------------------------------------------------------------
subroutine ConstructFillMatrix3DVectorial(PropertyID, EnterDataID, TimeID, &
HorizontalGridID, GeometryID, ExtractType, PointsToFill3D, &
Matrix3DU, Matrix3DV, Matrix3DW, Matrix3DX, Matrix3DY, &
TypeZUV, FillMatrix, FileNameHDF, ObjFillMatrix, &
OverrideValueKeyword, ClientID, PredictDTMethod, &
MinForDTDecrease, ValueIsUsedForDTPrediction, CheckDates, STAT)
!Arguments---------------------------------------------------------------
integer :: EnterDataID
integer :: TimeID
integer :: GeometryID
integer :: HorizontalGridID
integer :: ExtractType
integer, dimension(:, :, :), pointer :: PointsToFill3D
real, dimension(:, :, :), pointer :: Matrix3DU
real, dimension(:, :, :), pointer :: Matrix3DV
real, dimension(:, :, :), pointer :: Matrix3DX
real, dimension(:, :, :), pointer :: Matrix3DY
integer :: TypeZUV
type (T_PropertyID) :: PropertyID
real, dimension(:, :, :), pointer, optional :: Matrix3DW
integer, optional, intent(IN) :: ClientID
real , optional, intent(IN ) :: FillMatrix
character(*), dimension(2), optional, intent(IN ) :: FileNameHDF
character(*), optional, intent(IN ) :: OverrideValueKeyword
integer, optional, intent(INOUT) :: ObjFillMatrix
logical, optional, intent(IN) :: CheckDates
integer, optional, intent(OUT) :: STAT
integer, optional, intent(IN ) :: PredictDTMethod
real, optional, intent(IN ) :: MinForDTDecrease
logical, optional, intent(IN ) :: ValueIsUsedForDTPrediction
!Local-------------------------------------------------------------------
real :: FillMatrix_
integer :: ready_, STAT_, STAT_CALL, nUsers, ObjFillMatrix_
integer :: PredictDTMethod_
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mFillMatrix_)) then
nullify (FirstObjFillMatrix)
call RegisterModule (mFillMatrix_)
endif
if (present(ObjFillMatrix)) then
ObjFillMatrix_ = ObjFillMatrix
else
ObjFillMatrix_ = PropertyID%ObjFillMatrix
endif
call Ready(ObjFillMatrix_, ready_)
cd0 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
if (present(CheckDates)) then
Me%CheckDates = CheckDates
endif
Me%VectorialProp = .true.
!~ if (.not. Check_Vectorial_Property(PropertyID%IDNumber)) then
if (.not. PropertyID%IsVectorial) then
write(*,*) 'Constructing scalar property but expected vectorial'
stop 'ConstructFillMatrix3DVectorial - ModuleFillMatrix - ERR00'
endif
if (present(PredictDTMethod)) then
PredictDTMethod_ = PredictDTMethod
else
PredictDTMethod_ = 1
endif
if (present(MinForDTDecrease)) then
Me%MinForDTDecrease = MinForDTDecrease
endif
if (present(ValueIsUsedForDTPrediction)) then
Me%ValueIsUsedForDTPrediction = ValueIsUsedForDTPrediction
else
Me%ValueIsUsedForDTPrediction = .false.
endif
Me%ObjEnterData = AssociateInstance (mENTERDATA_, EnterDataID )
Me%ObjTime = AssociateInstance (mTIME_, TimeID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)
Me%ObjGeometry = AssociateInstance (mGEOMETRY_, GeometryID )
! JPW 2012-01-28: Use GeometrySize to set size of the matrix.
! Using ubound may cause inconsistency with padded matrices (making IUB larger than the actual grid bound).
! PointsToFill3D should have the same dimensions as Geometry anyway.
! Matrices that are used for the Thomas algorithm are padded if _USE_PAGELOCKED (CUDA only) or _PAD_MATRICES (Fortran) is defined
call GetGeometrySize(GeometryID , Me%Size3D, Me%WorkSize3D, STAT_)
! Me%Size3D%ILB = lbound(PointsToFill3D, dim = 1)
! Me%Size3D%IUB = ubound(PointsToFill3D, dim = 1)
! Me%Size3D%JLB = lbound(PointsToFill3D, dim = 2)
! Me%Size3D%JUB = ubound(PointsToFill3D, dim = 2)
! Me%Size3D%KLB = lbound(PointsToFill3D, dim = 3)
! Me%Size3D%KUB = ubound(PointsToFill3D, dim = 3)
! Me%WorkSize3D%ILB = Me%Size3D%ILB + 1
! Me%WorkSize3D%IUB = Me%Size3D%IUB - 1
! Me%WorkSize3D%JLB = Me%Size3D%JLB + 1
! Me%WorkSize3D%JUB = Me%Size3D%JUB - 1
! Me%WorkSize3D%KLB = Me%Size3D%KLB + 1
! Me%WorkSize3D%KUB = Me%Size3D%KUB - 1
Me%Size2D = T_Size2D(null_int, null_int, null_int, null_int)
Me%Dim = Dim3D
Me%TypeZUV = TypeZUV
Me%PropertyID = PropertyID
!Point to result matrixes
Me%Matrix3DU => Matrix3DU
Me%Matrix3DV => Matrix3DV
if(present(Matrix3DW)) Me%Matrix3DW => Matrix3DW
Me%PointsToFill3D => PointsToFill3D
!Point to input matrixes
Me%Matrix3DX => Matrix3DX
Me%Matrix3DY => Matrix3DY
!No need for Z that is the same as W (not transformation)
if (present(FillMatrix)) then
FillMatrix_ = FillMatrix
else
FillMatrix_ = null_real
endif
where (PointsToFill3D == WaterPoint) Me%Matrix3DU = FillMatrix_
where (PointsToFill3D == WaterPoint) Me%Matrix3DV = FillMatrix_
if(present(Matrix3DW)) where (PointsToFill3D == WaterPoint) Me%Matrix3DW = FillMatrix_
where (PointsToFill3D == WaterPoint) Me%Matrix3DX = FillMatrix_
where (PointsToFill3D == WaterPoint) Me%Matrix3DY = FillMatrix_
if (Me%TypeZUV == TypeU_) then
Me%Size3D%JUB = Me%Size3D%JUB + 1
Me%WorkSize3D%JUB = Me%WorkSize3D%JUB + 1
endif
if (Me%TypeZUV == TypeV_) then
Me%Size3D%IUB = Me%Size3D%IUB + 1
Me%WorkSize3D%IUB = Me%WorkSize3D%IUB + 1
endif
!Specific of the Vertical_Z matrix see ModuleGeometry
if (Me%TypeZUV == TypeW_) then
Me%WorkSize3D%KLB = Me%WorkSize3D%KLB - 1
endif
if (present(FileNameHDF)) then
Me%ArgumentFileName = .true.
Me%FileNameHDF(1) = trim(FileNameHDF(1))
Me%FileNameHDF(2) = trim(FileNameHDF(2))
else
Me%ArgumentFileName = .false.
endif
if(present(OverrideValueKeyword))then
Me%OverrideValueKeyword = trim(adjustl(OverrideValueKeyword))
Me%OverrideValueKeywordON = .true.
else
Me%OverrideValueKeywordON = .false.
end if
!Here matrixes Me%Matrix3DX and Me%Matrix3DY are filled from user info
if (present(ClientID)) then
call ReadOptions (ExtractType, &
PointsToFill3D = PointsToFill3D, &
PredictDTMethod = PredictDTMethod_, &
ClientID = ClientID)
else
call ReadOptions (ExtractType, &
PointsToFill3D = PointsToFill3D, &
PredictDTMethod = PredictDTMethod_)
endif
!!Need to rotate input field (Me%Matrix3DX and Me%Matrix3DY) to grid (Me%Matrix3DU and Me%Matrix3DV))
call RotateVectorFieldToGrid(HorizontalGridID = Me%ObjHorizontalGrid, &
VectorInX = Me%Matrix3DX, &
VectorInY = Me%Matrix3DY, &
VectorOutX = Me%Matrix3DU, &
VectorOutY = Me%Matrix3DV, &
WaterPoints3D = PointsToFill3D, &
RotateX = .true., &
RotateY = .true., &
KLB = Me%WorkSize3D%KLB, &
KUB = Me%WorkSize3D%KUB, &
STAT = STAT_CALL)
nUsers = DeassociateInstance(mENTERDATA_, Me%ObjEnterData)
if (nUsers == 0) stop 'ConstructFillMatrix3D - ModuleFillMatrix - ERR01'
if(Me%TimeEvolution == None)then
PropertyID%SolutionFromFile = .false.
else
PropertyID%SolutionFromFile = .true.
end if
!Returns ID
ObjFillMatrix_ = Me%InstanceID
if (present(ObjFillMatrix)) then
ObjFillMatrix = ObjFillMatrix_
else
PropertyID%ObjFillMatrix = ObjFillMatrix_
endif
!nullify will remove pointer to Matrix3DU and Matrix3DV (Module field) since they were alrady modifed
!when modifying Me%Matrix3DU and Me%Matrix3DV
nullify(Me%Matrix3D )
nullify(Me%Matrix3DU )
nullify(Me%Matrix3DV )
nullify(Me%Matrix3DW )
nullify(Me%Matrix3DX )
nullify(Me%Matrix3DY )
nullify(Me%PointsToFill3D)
STAT_ = SUCCESS_
else cd0
stop 'ConstructFillMatrix2D - ModuleFillMatrix - ERR02'
end if cd0
if (present(STAT)) STAT = STAT_
!----------------------------------------------------------------------
end subroutine ConstructFillMatrix3DVectorial
!----------------------------------------------------------------------
subroutine AllocateInstance
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
type (T_FillMatrix), pointer :: NewObjFillMatrix
type (T_FillMatrix), pointer :: PreviousObjFillMatrix
!Allocates new instance
allocate (NewObjFillMatrix)
nullify (NewObjFillMatrix%Next)
!Insert New Instance into list and makes Current point to it
if (.not. associated(FirstObjFillMatrix)) then
FirstObjFillMatrix => NewObjFillMatrix
Me => NewObjFillMatrix
else
PreviousObjFillMatrix => FirstObjFillMatrix
Me => FirstObjFillMatrix%Next
do while (associated(Me))
PreviousObjFillMatrix => Me
Me => Me%Next
enddo
Me => NewObjFillMatrix
PreviousObjFillMatrix%Next => NewObjFillMatrix
endif
Me%InstanceID = RegisterNewInstance (mFILLMATRIX_)
end subroutine AllocateInstance
!--------------------------------------------------------------------------
subroutine ReadOptions (ExtractType, PointsToFill2D, PointsToFill3D, ClientID, PredictDTMethod)
!Arguments-------------------------------------------------------------
integer :: ExtractType
integer, dimension(:, :), pointer, optional :: PointsToFill2D
integer, dimension(:, :, :), pointer, optional :: PointsToFill3D
integer, optional, intent(IN) :: ClientID
integer, optional, intent(IN) :: PredictDTMethod
!Local----------------------------------------------------------------
integer :: STAT_CALL, ClientID_
integer :: iflag
character(len=StringLength) :: AuxString
logical :: AuxBoolean
integer :: AuxInteger
logical :: CanContinue
!---------------------------------------------------------------------
if (present(ClientID)) then
ClientID_ = ClientID
else
ClientID_ = FillValueInt
endif
!Reads Time Evolution
call GetData(AuxString, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'FILE_IN_TIME', &
default = "None", &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR010'
!write(*,*) trim(adjustl(AuxString))
select case (trim(adjustl(AuxString)))
case ("None", "NONE", "none")
Me%TimeEvolution = None
case ("Hdf", "HDF", "hdf")
Me%TimeEvolution = ReadHDF
case ("Timeserie", "TIMESERIE", "timeserie", "TimeSerie")
Me%TimeEvolution = ReadTimeSerie
case ("Profile_Timeserie", "PROFILE_TIMESERIE", "profile_timeserie", "Profile_TimeSerie")
Me%TimeEvolution = ProfileTimeSerie
case ("Multitimeserie", "MULTITIMESERIE", "multitimeserie", "MultiTimeserie")
Me%TimeEvolution = MultiTimeserie
case ("Analytic Wave", "ANALYTIC WAVE", "analytic wave")
Me%TimeEvolution = AnalyticWave
case default
write(*,*)'Invalid option for keyword FILE_IN_TIME'
stop 'ReadOptions - ModuleFillMatrix - ERR020'
end select
!write (*,*) Me%TimeEvolution
if (Me%ArgumentFileName) Me%TimeEvolution = ReadHDF
!write (*,*) Me%TimeEvolution
if(Me%TimeEvolution == None)then
call GetData(AuxString, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'INITIALIZATION_METHOD', &
default = "Constant", &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR030'
!write(*,*) trim(adjustl(AuxString))
select case (trim(adjustl(AuxString)))
case ("Constant", "CONSTANT", "constant")
Me%InitializationMethod = Constant
case ("Layers", "LAYERS", "layers")
Me%InitializationMethod = Layers
case ("Boxes", "BOXES", "boxes")
Me%InitializationMethod = Boxes
case ("ASCII_File", "ASCII_FILE", "ascii_file", "Ascii_file")
Me%InitializationMethod = AsciiFile
case ("Profile", "PROFILE", "profile")
Me%InitializationMethod = Profile
case ("Analytic Profile", "ANALYTIC PROFILE", "analytic profile")
Me%InitializationMethod = AnalyticProfile
case ("Hdf", "HDF", "hdf")
Me%InitializationMethod = ReadHDF
case ("Timeserie", "TIMESERIE", "timeserie", "TimeSerie")
Me%InitializationMethod = ReadTimeSerie
case ("Profile_Timeserie", "PROFILE_TIMESERIE", "profile_timeserie", "Profile_TimeSerie")
Me%InitializationMethod = ProfileTimeSerie
case ("Sponge", "SPONGE", "sponge")
Me%InitializationMethod = Sponge
case default
write(*,*)'Invalid option for keyword INITIALIZATION_METHOD'
stop 'ReadOptions - ModuleFillMatrix - ERR040'
end select
call GetData(Me%RemainsConstant, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'REMAIN_CONSTANT', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR050'
call GetData(AuxString, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'INITIALIZATION_DEFAULT', &
default = "Constant", &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR060'
select case (trim(adjustl(AuxString)))
case ("Constant", "CONSTANT", "constant")
Me%InitializationDefault = Constant
case ("Profile_Timeserie", "PROFILE_TIMESERIE", "profile_timeserie", "Profile_TimeSerie")
Me%InitializationDefault = ProfileTimeSerie
end select
end if
! Check if the simulation goes backward in time or forward in time (default mode)
call GetBackTracking(Me%ObjTime, Me%BackTracking, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR070'
call GetComputeTimeLimits(Me%ObjTime, BeginTime = Me%BeginTime, &
EndTime = Me%EndTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR080'
!Gets the default value
call GetData(Me%DefaultValue, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'DEFAULTVALUE', &
ClientModule = 'ModuleFillMatrix', &
default = null_real, &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) then
!normal case, found only one value (STAT_CALL size err and iflag = 1). ignore
!vectorial 2D case, found only two values (STAT_CALL size err and iflag = 2). ignore
if (STAT_CALL /= SIZE_ERR_) then
write(*,*) "Property ", trim(Me%PropertyID%Name)
write(*,*) "Wrong number of arguments for DEFAULTVALUE keyword"
stop 'ReadOptions - ModuleFillMatrix - ERR090'
endif
!need to define two or three values, not only one
if (Me%VectorialProp .and. iflag == 1) then
write(*,*) 'DEFAULTVALUE for vectorial prop needs to define the x and y'
write(*,*) '(and optionally z) components in 2/3 values'
write(*,*) 'separated by space. Property : ', trim(Me%PropertyID%Name)
stop 'ReadOptions - ModuleFillMatrix - ERR095'
endif
endif
if (iflag == 0 .and. .not. Me%ArgumentFileName) then
if(Me%OverrideValueKeywordON .and. Me%InitializationMethod == Constant)then
call GetData(Me%DefaultValue, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = trim(Me%OverrideValueKeyword), &
ClientModule = 'ModuleFillMatrix', &
default = null_real, &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) then
if(STAT_CALL == SIZE_ERR_)then
if(.not. Me%VectorialProp .and. iflag .ne. 1)then
write(*,*) "Property ", trim(Me%PropertyID%Name)
write(*,*) "DEFAULTVALUE must be a scalar"
stop 'ReadOptions - ModuleFillMatrix - ERR99'
elseif(iflag > 3)then
write(*,*) "Property ", trim(Me%PropertyID%Name)
write(*,*) "Too many arguments in keyword "//trim(Me%OverrideValueKeyword)
stop 'ReadOptions - ModuleFillMatrix - ERR98'
endif
else
stop 'ReadOptions - ModuleFillMatrix - ERR100'
endif
end if
if(iflag == 0)then
write(*,*)'Please define override keyword '//trim(Me%OverrideValueKeyword)
write(*,*)'to give a default value for property '//trim(Me%PropertyID%Name)
stop 'ReadOptions - ModuleFillMatrix - ERR101'
end if
elseif(Me%OverrideValueKeywordON .and. Me%InitializationMethod == Boxes)then
Me%DefaultValue = null_real
elseif(Me%OverrideValueKeywordON .and. Me%InitializationMethod == AsciiFile)then
Me%DefaultValue = null_real
elseif((Me%OverrideValueKeywordON .and. .not. Me%InitializationMethod == Boxes ) .or. &
(Me%OverrideValueKeywordON .and. .not. Me%InitializationMethod == AsciiFile) .or. &
(Me%OverrideValueKeywordON .and. .not. Me%InitializationMethod == Constant ))then
write(*,*)'Initialization method for property '//trim(Me%PropertyID%Name)
write(*,*)'can only be CONSTANT, BOXES or ASCII'
stop 'ReadOptions - ModuleFillMatrix - ERR102'
else
write(*,*)'Please define default value for property '//trim(Me%PropertyID%Name)
stop 'ReadOptions - ModuleFillMatrix - ERR103'
end if
end if
!Keywords that define how to "work" values (interpolate, accumulate, nothing)
call GetData(Me%InterpolateValues, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'INTERPOLATE_VALUES', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR104'
call GetData(Me%AccumulateValues, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'ACCUMULATE_VALUES', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR105'
call GetData(Me%UseOriginalValues, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'USE_ORIGINAL_VALUES', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR110'
if ((Me%InterpolateValues .AND. Me%AccumulateValues) .OR. &
(Me%InterpolateValues .AND. Me%UseOriginalValues) .OR. &
(Me%AccumulateValues .AND. Me%UseOriginalValues)) then
write (*,*) 'The keywords INTERPOLATE_VALUES, ACCUMULATE_VALUES and'
write (*,*) 'USE_ORIGINAL_VALUES are mutually exclusives.'
write (*,*) 'Only one can be set to true.'
stop 'ReadOptions - ModuleFillMatrix - ERR120'
endif
!Property ValuesType (Interpolated, accumulate, original value)
call GetData(AuxInteger, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'VALUES_TYPE', &
default = InterpolatedValues, &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR120-A'
if (iflag .NE. 0) then
write(*,*)
write(*,*) 'ModuleFillMatrix WARNING:'
write(*,*) 'VALUES_TYPE keyword was removed.'
write(*,*) 'Use these instead: '
write(*,*) ' INTERPOLATE_VALUES => to interpolate values'
write(*,*) ' ACCUMULATE_VALUES => to accumulate values'
write(*,*) ' USE_ORIGINAL_VALUES => to use original values'
write(*,*)
stop 'ReadOptions - ModuleFillMatrix - ERR120-B'
endif
!Property not interpolated (e.g. Concentration on rain)
call GetData(AuxBoolean, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'NO_INTERPOLATION_OR_ACCUMULATION', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR130'
if (iflag .NE. 0) then
write(*,*)
write(*,*) 'ModuleFillMatrix WARNING:'
write(*,*) 'NO_INTERPOLATION_OR_ACCUMULATION keyword is deprecated.'
write(*,*) 'To use values without interpolation or accumulation use instead: '
write(*,*) ' "USE_ORIGINAL_VALUES : 1"'
write(*,*)
if (AuxBoolean) then
Me%InterpolateValues = .false.
Me%AccumulateValues = .false.
Me%UseOriginalValues = .true.
endif
endif
if (.NOT. AuxBoolean) then
!Accumulitve Property (e.g. Rain from gauges) !this keyword name should be changed because is
!not describing well what is the computation. However a lot of people use it already...
call GetData(AuxBoolean, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'NO_INTERPOLATION', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR140'
if (iflag .NE. 0) then
write(*,*)
write(*,*) 'ModuleFillMatrix WARNING:'
write(*,*) 'NO_INTERPOLATION keyword is deprecated.'
write(*,*) 'To use interpolated values use instead : '
write(*,*) ' "INTERPOLATE_VALUES : 1"'
write(*,*) 'To use accumulated values use instead :'
write(*,*) ' "ACCUMULATE_VALUES : 1" '
write(*,*)
if (AuxBoolean) then
Me%InterpolateValues = .false.
Me%AccumulateValues = .true.
Me%UseOriginalValues = .false.
else
Me%InterpolateValues = .true.
Me%AccumulateValues = .false.
Me%UseOriginalValues = .false.
endif
endif
endif
if ((.NOT. Me%InterpolateValues) .AND. &
(.NOT. Me%AccumulateValues) .AND. &
(.NOT. Me%UseOriginalValues)) then
Me%InterpolateValues = .true.
endif
!When to shut down DT?
call GetData(Me%MinForDTDecrease, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'MIN_FOR_DT_DECREASE', &
default = AllMostZero, &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR150'
!Assumes the last instant value to avoid linear interpolation
call GetData(Me%PreviousInstantValues, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'PREVIOUS_INSTANT_VALUES', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR155'
call GetData(Me%IgnoreNoDataPoint, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'IGNORE_NODATA_POINT', &
default = .true., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR156'
if (.NOT. Me%IgnoreNoDataPoint) then
call GetData(Me%NoDataValue, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'NODATA_VALUE', &
default = -99.0, &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR157'
endif
call GetData(Me%PredictDTMethod, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'PREDICT_DT_METHOD', &
default = PredictDTMethod, &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ReadOptions - ModuleFillMatrix - ERR158'
!Fill Matrix with default Value. This is the user field (not rotated to grid)
if (Me%VectorialProp) then
if (Me%Dim == Dim2D) then
where (PointsToFill2D == WaterPoint) Me%Matrix2DX = Me%DefaultValue(1)
where (PointsToFill2D == WaterPoint) Me%Matrix2DY = Me%DefaultValue(2)
else
where (PointsToFill3D == WaterPoint) Me%Matrix3DX = Me%DefaultValue(1)
where (PointsToFill3D == WaterPoint) Me%Matrix3DY = Me%DefaultValue(2)
endif
!User field angle (not rotated to grid)
else if (Me%RotateAngleToGrid) then
if (Me%Dim == Dim2D) then
where (PointsToFill2D == WaterPoint) Me%Matrix2DFieldAngle = Me%DefaultValue(1)
else
where (PointsToFill3D == WaterPoint) Me%Matrix3DFieldAngle = Me%DefaultValue(1)
endif
else
if (Me%Dim == Dim2D) then
where (PointsToFill2D == WaterPoint) Me%Matrix2D = Me%DefaultValue(1)
else
where (PointsToFill3D == WaterPoint) Me%Matrix3D = Me%DefaultValue(1)
endif
endif
CanContinue = .true.
if (Me%VectorialProp) then
if (Me%TimeEvolution == None) then
if ((Me%InitializationDefault /= Constant) .or. (Me%InitializationMethod /= Constant &
.and. Me%InitializationMethod /= ASCIIFile .and. Me%InitializationMethod /= ReadTimeSerie &
.and. Me%InitializationMethod /= ReadHDF)) then
CanContinue = .false.
endif
elseif (Me%TimeEvolution /= ReadTimeSerie .and. Me%TimeEvolution /= ReadHDF) then
CanContinue = .false.
endif
endif
if (.not. CanContinue) then
write(*,*) 'For now vectorial property can only be defined by constant, asciifile, timeserie or HDF'
write(*,*) 'Property : ', trim(Me%PropertyID%Name)
stop ' ReadOptios - ModuleFillMatrix - ERR159'
endif
select case (Me%TimeEvolution)
case(None)
!Fill Matrix with an initialization default
select case(Me%InitializationDefault)
case(Constant)
!Do nothing. Matrix values equal default value
case(ProfileTimeSerie)
if (Me%Dim == Dim2D) then
write(*,*)'Cannot Initialize 2D matrix by profile time serie'
stop 'ReadOptions - ModuleFillMatrix - ERR160'
else
call ConstructProfileTSDefault (PointsToFill3D, ExtractType)
endif
end select
select case(Me%InitializationMethod)
case(Constant)
!Do nothing. Matrix values equal default value
case(Layers)
if (Me%Dim == Dim2D) then
write(*,*)'Cannot Initialise 2D Matrix by Layers'
stop 'ReadOptions - ModuleFillMatrix - ERR170'
else
call ConstructSpaceLayers (ExtractType, PointsToFill3D)
endif
case(Boxes)
if (Me%Dim == Dim2D) then
call ConstructSpaceBox (ExtractType, PointsToFill2D = PointsToFill2D)
else
call ConstructSpaceBox (ExtractType, PointsToFill3D = PointsToFill3D)
endif
case(AsciiFile)
if (Me%Dim == Dim2D) then
call ConstructSpaceASCIIFile (ExtractType, PointsToFill2D = PointsToFill2D)
else
call ConstructSpaceASCIIFile (ExtractType, PointsToFill3D = PointsToFill3D)
endif
case(Sponge)
if (Me%Dim == Dim2D) then
call ConstructSponge (ExtractType, PointsToFill2D = PointsToFill2D)
else
call ConstructSponge (ExtractType, PointsToFill3D = PointsToFill3D)
endif
case(Profile)
if (Me%Dim == Dim2D) then
write(*,*)'Cannot Initialise 2D Matrix by Profile'
stop 'ReadOptions - ModuleFillMatrix - ERR180'
else
call ConstructSpaceProfile(ExtractType, PointsToFill3D)
endif
case(AnalyticProfile)
if (Me%Dim == Dim2D) then
write(*,*)'Cannot Initialise 2D Matrix by Profile'
stop 'ReadOptions - ModuleFillMatrix - ERR190'
else
call ConstructAnalyticProfile(ExtractType, PointsToFill3D)
endif
case(ReadTimeSerie)
if (Me%Dim == Dim2D) then
call ConstructSpaceTimeSerie (ExtractType, PointsToFill2D = PointsToFill2D)
else
call ConstructSpaceTimeSerie (ExtractType, PointsToFill3D = PointsToFill3D)
endif
case(ReadHDF)
if (Me%Dim == Dim2D) then
call ConstructHDFInput (ExtractType, ClientID_, PointsToFill2D = PointsToFill2D)
else
call ConstructHDFInput (ExtractType, ClientID_, PointsToFill3D = PointsToFill3D)
endif
case(ProfileTimeSerie)
call ConstructProfileTimeSerie (PointsToFill3D, ExtractType)
end select
case(ReadTimeSerie)
if (Me%Dim == Dim2D) then
call ConstructSpaceTimeSerie (ExtractType, PointsToFill2D = PointsToFill2D)
else
call ConstructSpaceTimeSerie (ExtractType, PointsToFill3D = PointsToFill3D)
endif
case(ReadHDF)
if (Me%Dim == Dim2D) then
call ConstructHDFInput (ExtractType, ClientID_, PointsToFill2D = PointsToFill2D)
else
call ConstructHDFInput (ExtractType, ClientID_, PointsToFill3D = PointsToFill3D)
endif
case(ProfileTimeSerie)
if (Me%Dim == Dim2D) then
write(*,*)'Cannot read 2D matrix by profile time serie'
stop 'ReadOptions - ModuleFillMatrix - ERR200'
else
call ConstructProfileTimeSerie (PointsToFill3D, ExtractType)
endif
case(MultiTimeserie)
if (.not. present(ClientID)) &
stop 'ReadOptions - ModuleFillMatrix - ERR210'
if (Me%Dim == Dim2D) then
call ConstructMultiTimeSerie (ClientID, PointsToFill2D = PointsToFill2D)
call ModifyMultiTimeSerie (PointsToFill2D = PointsToFill2D)
else
call ConstructMultiTimeSerie (ClientID, PointsToFill3D = PointsToFill3D)
call ModifyMultiTimeSerie (PointsToFill3D = PointsToFill3D)
endif
case(AnalyticWave)
call ConstructAnalyticWave (ExtractType)
if (Me%Dim == Dim2D) then
call ModifyAnalyticWave (PointsToFill2D = PointsToFill2D)
else
call ModifyAnalyticWave (PointsToFill3D = PointsToFill3D)
endif
end select
end subroutine ReadOptions
!--------------------------------------------------------------------------
subroutine InsertTimeSerieToList(TimeSerie)
!Arguments-------------------------------------------------------------
type (T_TimeSerie), pointer :: TimeSerie
!Local-----------------------------------------------------------------
type (T_TimeSerie), pointer :: CurrentTimeSerie => null()
type (T_TimeSerie), pointer :: PreviousTimeSerie => null()
!Inserts a new property to the list of properties
if (.not. associated(Me%FirstTimeSerie)) then
Me%FirstTimeSerie => TimeSerie
CurrentTimeSerie => TimeSerie
else
PreviousTimeSerie => Me%FirstTimeSerie
CurrentTimeSerie => PreviousTimeSerie%Next
do while (associated(CurrentTimeSerie))
PreviousTimeSerie => CurrentTimeSerie
CurrentTimeSerie => PreviousTimeSerie%Next
enddo
PreviousTimeSerie%Next => TimeSerie
TimeSerie%Prev => PreviousTimeSerie
endif
Me%nTimeSeries = Me%nTimeSeries + 1
end subroutine InsertTimeSerieToList
!--------------------------------------------------------------------------
subroutine InsertASCIIFileToList(ASCIIFile)
!Arguments-------------------------------------------------------------
type (T_ASCIIFile), pointer :: ASCIIFile
!Local-----------------------------------------------------------------
type (T_ASCIIFile), pointer :: CurrentASCIIFile => null()
type (T_ASCIIFile), pointer :: PreviousASCIIFile => null()
!Inserts a new property to the list of properties
if (.not. associated(Me%FirstASCIIFile)) then
Me%FirstASCIIFile => ASCIIFile
CurrentASCIIFile => ASCIIFile
else
PreviousASCIIFile => Me%FirstASCIIFile
CurrentASCIIFile => PreviousASCIIFile%Next
do while (associated(CurrentASCIIFile))
PreviousASCIIFile => CurrentASCIIFile
CurrentASCIIFile => PreviousASCIIFile%Next
enddo
PreviousASCIIFile%Next => ASCIIFile
ASCIIFile%Prev => PreviousASCIIFile
endif
Me%nASCIIFiles = Me%nASCIIFiles + 1
end subroutine InsertASCIIFileToList
!--------------------------------------------------------------------------
subroutine InsertHDFToList(HDF)
!Arguments-------------------------------------------------------------
type (T_Field4D), pointer :: HDF
!Local-----------------------------------------------------------------
type (T_Field4D), pointer :: CurrentHDF => null()
type (T_Field4D), pointer :: PreviousHDF => null()
!Inserts a new property to the list of properties
if (.not. associated(Me%FirstHDF)) then
Me%FirstHDF => HDF
CurrentHDF => HDF
else
PreviousHDF => Me%FirstHDF
CurrentHDF => PreviousHDF%Next
do while (associated(CurrentHDF))
PreviousHDF => CurrentHDF
CurrentHDF => PreviousHDF%Next
enddo
PreviousHDF%Next => HDF
HDF%Prev => PreviousHDF
endif
Me%nHDFs = Me%nHDFs + 1
end subroutine InsertHDFToList
!--------------------------------------------------------------------------
subroutine ConstructMultiTimeSerie (ClientID, PointsToFill2D, PointsToFill3D)
!Arguments-------------------------------------------------------------
integer, intent(IN) :: ClientID
integer, dimension(:, :), pointer, optional :: PointsToFill2D
integer, dimension(:, :, :), pointer, optional :: PointsToFill3D
!Local----------------------------------------------------------------
character(PathLength) :: FileName
integer :: MaskGridDataID = 0
integer :: TypeZUV
real, dimension(:, :), pointer :: GridData2D
real, dimension(:, :, :), pointer :: GridData3D
integer :: STAT_CALL
integer :: iflag
integer :: ilb, iub, jlb, jub, klb, kub
integer :: i, j, k
integer :: index
logical :: found
type(T_Station), dimension(:), pointer :: sl
integer :: ClientNumber
Type(T_Time) :: CurrentTime
!----------------------------------------------------------------------
ClientNumber = ClientID
!Gets the number of Source blocks
call GetNumberOfBlocks(Me%ObjEnterData, BeginMTSBlock, EndMTSBlock, &
FromBlock_, Me%MultiTimeSerie%NumberOfSources, &
ClientNumber, STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR010'
allocate (Me%MultiTimeSerie%StationsList(Me%MultiTimeSerie%NumberOfSources))
sl => Me%MultiTimeSerie%StationsList
do index = 1, Me%MultiTimeSerie%NumberOfSources
call ExtractBlockFromBlock(Me%ObjEnterData, &
ClientNumber = ClientNumber, &
block_begin = BeginMTSBlock, &
block_end = EndMTSBlock, &
BlockInBlockFound = found, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR020'
if (found) then
call GetData(sl(index)%FillID, &
Me%ObjEnterData , iflag, &
SearchType = FromBlockInBlock, &
keyword = 'FILL_ID', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR021'
if (iflag /= 1) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR022'
call GetData(sl(index)%RemainConstant, &
Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlock, &
keyword = 'REMAIN_CONSTANT', &
default = .false., &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR023'
if (sl(index)%RemainConstant) then
call GetData(sl(index)%NewValue, &
Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlock, &
keyword = 'DEFAULTVALUE', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR024'
else
call GetData(sl(index)%FileName, &
Me%ObjEnterData, iflag, &
SearchType = FromBlockInBlock, &
keyword = 'FILE_NAME', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR025'
call GetData(sl(index)%Column, &
Me%ObjEnterData , iflag, &
SearchType = FromBlockInBlock, &
keyword = 'DATA_COLUMN', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR026'
if (iflag /= 1) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR027'
!Starts Time Serie
call StartTimeSerieInput(sl(index)%ObjTimeSerie, &
sl(index)%FileName, &
Me%ObjTime, &
CheckDates = Me%CheckDates, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR028'
endif
else
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR030'
endif
enddo
!call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT_CALL)
!if (STAT_CALL /= SUCCESS_) &
! stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR040'
!Gets the name of the mask file
call GetData(FileName, &
Me%ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'ASCII_MASK_FILENAME', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR050'
if (iflag /= 1)then
write(*,*)'ASCII_MASK_FILENAME is missing'
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR051'
end if
call GetData(Me%MultiTimeSerie%DataProcessing, &
Me%ObjEnterData , iflag, &
SearchType = FromBlock, &
keyword = 'DATA_PROCESSING', &
Default = Interpolate, &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR060'
if (iflag /= 1) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR070'
if (Me%Dim == Dim2D) then
ilb = Me%WorkSize2D%ILB
iub = Me%WorkSize2D%IUB
jlb = Me%WorkSize2D%JLB
jub = Me%WorkSize2D%JUB
allocate (Me%MultiTimeSerie%FillGrid2D(ilb:iub, jlb:jub))
call ConstructGridData(MaskGridDataID, Me%ObjHorizontalGrid, &
FileName = FileName, &
DefaultValue = Me%DefaultValue(1), &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR080'
call GetGridDataType(MaskGridDataID, TypeZUV, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR090'
if(TypeZUV .ne. Me%TypeZUV)then
write(*,*)'Inconsistency found in type ZUV'
write(*,*)'Grid data: '//trim(FileName)
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR100'
end if
call GetGridData (MaskGridDataID, GridData2D, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR110'
do j = jlb, jub
do i = ilb, iub
if (PointsToFill2D(i,j) == WaterPoint) then
found = .false.
do index = 1, Me%MultiTimeSerie%NumberOfSources
if (sl(index)%FillID == GridData2D(i,j)) then
found = .true.
exit
endif
enddo
if (found) then
Me%MultiTimeSerie%FillGrid2D(i,j) = index
else
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR120'
endif
else
Me%MultiTimeSerie%FillGrid2D(i,j) = -99
endif
enddo
enddo
call UnGetGridData (MaskGridDataID, GridData2D, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR130'
else
ilb = Me%WorkSize3D%ILB
iub = Me%WorkSize3D%IUB
jlb = Me%WorkSize3D%JLB
jub = Me%WorkSize3D%JUB
klb = Me%WorkSize3D%KLB
kub = Me%WorkSize3D%KUB
allocate (Me%MultiTimeSerie%FillGrid3D(ilb:iub, jlb:jub, klb:kub))
call ConstructGridData(MaskGridDataID, Me%ObjHorizontalGrid, &
FileName = FileName, &
KLB = Me%Worksize3D%KLB, &
KUB = Me%Worksize3D%KUB, &
DefaultValue = Me%DefaultValue(1), &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR140'
call GetGridDataType(Me%ASCIIFile%GridDataID, TypeZUV, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR150'
if(TypeZUV .ne. Me%TypeZUV)then
write(*,*)'Inconsistency found in type ZUV'
write(*,*)'Grid data: '//trim(FileName)
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR160'
end if
call GetGridData (MaskGridDataID, GridData3D, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR170'
do j = jlb, jub
do i = ilb, iub
do k = klb, kub
if (PointsToFill3D(i,j,k) == WaterPoint) then
found = .false.
do index = 1, Me%MultiTimeSerie%NumberOfSources
if (sl(index)%FillID == GridData3D(i,j,k)) then
found = .true.
exit
endif
enddo
if (found) then
Me%MultiTimeSerie%FillGrid3D(i,j,k) = index
else
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR180'
endif
else
Me%MultiTimeSerie%FillGrid3D(i,j,k) = -99
endif
enddo
enddo
enddo
call UnGetGridData (MaskGridDataID, GridData3D, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR190'
endif
call KillGridData (MaskGridDataID, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR200'
if (Me%PredictDTMethod == 2) then
!Gets Current Time
call GetComputeCurrentTime(Me%ObjTime, CurrentTime, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR210'
Me%DTForNextEvent = -null_real
Me%PredictedDT = -null_real
Me%DTForNextDataset = -null_real
Me%NextValueForDTPred = 0.0
sl => Me%MultiTimeSerie%StationsList
do index = 1, Me%MultiTimeSerie%NumberOfSources
if (.not. sl(index)%RemainConstant) then
call GetTimeSerieDataValues(sl(index)%ObjTimeSerie, &
sl(index)%NumberOfInstants, &
STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR219'
if(sl(index)%NumberOfInstants > 1)then
sl(index)%PreviousInstant = 1
call GetTimeSerieTimeOfDataset(sl(index)%ObjTimeSerie, &
sl(index)%PreviousInstant, &
sl(index)%PreviousTime, &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR220'
sl(index)%NextInstant = 1
sl(index)%NextTime = sl(index)%PreviousTime
call ActualizeMultiTimeSerieValues (sl(index))
sl(index)%NextEventStart = sl(index)%PreviousTime
sl(index)%NextEventEnd = sl(index)%PreviousTime
sl(index)%NextValueForDTPred = sl(index)%NextValue
sl(index)%DTForNextEvent = 0.0
!
! if (Me%AccumulateValues) then
! if (sl(index)%NextValue > Me%MinForDTDecrease) then
! sl(index)%NextEventStart = sl(index)%PreviousTime
! sl(index)%NextEventEnd = sl(index)%NextTime
! else
! call FindNextEventInMultiTimeSerie(CurrentTime, sl(index))
! endif
! else
! sl(index)%DTForNextEvent = 0.0
! endif
!
! if (Me%ValueIsUsedForDTPrediction .and. (sl(index)%DTForNextEvent == 0.0)) then
! if (Me%UseOriginalValues) then
! if (sl(index)%NextValue > sl(index)%NextValueForDTPred) &
! sl(index)%NextValueForDTPred = sl(index)%NextValue
! elseif (Me%AccumulateValues) then
! value = sl(index)%NextValue / (sl(index)%PreviousTime - sl(index)%NextTime)
! if (value > Me%NextValueForDTPred) &
! sl(index)%NextValueForDTPred = value
! else
! stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR240'
! endif
! endif
else
stop 'ConstructMultiTimeSerie - ModuleFillMatrix - ERR241'
endif
else
if (sl(index)%NewValue > 0.0) then
sl(index)%DTForNextEvent = 0.0
sl(index)%PredictedDT = -null_real
sl(index)%DTForNextDataset = -null_real
sl(index)%NextValueForDTPred = sl(index)%NewValue
else
sl(index)%DTForNextEvent = -null_real
sl(index)%PredictedDT = -null_real
sl(index)%DTForNextDataset = -null_real
sl(index)%NextValueForDTPred = 0.0
endif
endif
!
! if (Me%DTForNextEvent > sl(index)%DTForNextEvent) &
! Me%DTForNextEvent = sl(index)%DTForNextEvent
!
! if (Me%PredictedDT > sl(index)%PredictedDT) &
! Me%PredictedDT = sl(index)%PredictedDT
!
! if (Me%DTForNextDataset > sl(index)%DTForNextDataset) &
! Me%DTForNextDataset = sl(index)%DTForNextDataset
!
! if ((sl(index)%DTForNextEvent == 0.0) .and. &
! (Me%NextValueForDTPred < sl(index)%NextValueForDTPred)) then
! Me%NextValueForDTPred = sl(index)%NextValueForDTPred
! endif
enddo
endif
!----------------------------------------------------------------------
end subroutine ConstructMultiTimeSerie
!--------------------------------------------------------------------------
subroutine ConstructSpaceLayers (ExtractType, PointsToFill3D)
!Arguments-------------------------------------------------------------
integer :: ExtractType
integer, dimension(:, :, :), pointer :: PointsToFill3D
!Local----------------------------------------------------------------
character(PathLength) :: FileName
integer :: STAT_CALL
integer :: iflag
integer :: k, NLayers
integer :: ObjEnterData = 0
integer :: ClientNumber
integer :: FirstLine, LastLine
logical :: BlockFound
integer :: line, l
real, dimension(:), allocatable :: Aux
!----------------------------------------------------------------------
!Begin----------------------------------------------------------------
allocate (Me%Layers%Values(Me%WorkSize3D%KUB))
call GetData(FileName, &
Me%ObjEnterData , iflag, &
SearchType = ExtractType, &
keyword = 'FILENAME', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR01'
flag0: if (iflag /= 0) then
!Opens File
call ConstructEnterData(ObjEnterData, FileName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR02'
call ExtractBlockFromBuffer(ObjEnterData, ClientNumber, &
BeginLayers, EndLayers, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR03'
BF: if (BlockFound) then
NLayers = LastLine - FirstLine - 1
if (NLayers > Me%WorkSize3D%KUB) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR04'
!Allocates auxiliar variables
allocate (Aux(2))
l = 1
do line = FirstLine + 1, LastLine - 1
call GetData(Aux, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = line, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR05'
if (Aux(1) < Me%WorkSize3D%KLB .or. Aux(1) > Me%WorkSize3D%KUB) &
stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR06'
Me%Layers%Values(int(Aux(1))) = Aux(2)
l = l + 1
enddo
deallocate(Aux)
call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR07'
call KillEnterData(ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR08'
else
write(*,*) 'Block <BeginLayers>, <EndLayers> not found'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR09'
endif BF
else if (iflag == 0) then flag0
call GetData(Me%Layers%Values, &
Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'LAYERS_VALUES', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR10'
if (iflag /= Me%WorkSize3D%KUB) then
write(*,*)'Invalid Number of Layers'
stop 'ConstructSpaceLayers - ModuleFillMatrix - ERR11'
endif
endif flag0
!Fill Matrix with default Value
do k = Me%WorkSize3D%KLB, Me%WorkSize3D%KUB
where (PointsToFill3D(:,:,k) == WaterPoint) Me%Matrix3D(:,:,k) = Me%Layers%Values(k)
enddo
end subroutine ConstructSpaceLayers
!--------------------------------------------------------------------------
subroutine ConstructSpaceProfile (ExtractType, PointsToFill3D)
!Arguments-------------------------------------------------------------
integer :: ExtractType
integer, dimension(:, :, :), pointer :: PointsToFill3D
!Local----------------------------------------------------------------
real, dimension(:,:,:), pointer :: SZZ
real, dimension(:), pointer :: Values, Depth
real :: CellDepth
character(PathLength) :: FileName
integer :: STAT_CALL
integer :: iflag
integer :: i, j, k, NDEPTHS
integer :: ILB, IUB, JLB, JUB, KLB, KUB
integer :: ObjEnterData = 0
integer :: ClientNumber
integer :: FirstLine, LastLine
logical :: BlockFound
integer :: line, l
real, dimension(:), allocatable :: Aux
!----------------------------------------------------------------------
!Begin----------------------------------------------------------------
ILB = Me%WorkSize3D%ILB
IUB = Me%WorkSize3D%IUB
JLB = Me%WorkSize3D%JLB
JUB = Me%WorkSize3D%JUB
KLB = Me%WorkSize3D%KLB
KUB = Me%WorkSize3D%KUB
!Gets Center of the cells
call GetGeometryDistances(Me%ObjGeometry, SZZ = SZZ, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR01'
call GetData(FileName, &
Me%ObjEnterData , iflag, &
SearchType = ExtractType, &
keyword = 'FILENAME', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR02'
flag0: if (iflag /= 0) then
!Opens File
call ConstructEnterData(ObjEnterData, FileName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR03'
call ExtractBlockFromBuffer(ObjEnterData, ClientNumber, &
BeginProfile, EndProfile, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR04'
BF: if (BlockFound) then
NDEPTHS = LastLine - FirstLine - 1
!Allocates auxiliar variables
allocate (Values (NDEPTHS))
allocate (Depth (NDEPTHS))
allocate (Aux(2))
l = 1
do line = FirstLine + 1, LastLine - 1
call GetData(Aux, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = line, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR05'
Depth (l) = Aux(1)
Values(l) = Aux(2)
l = l + 1
enddo
deallocate(Aux)
call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR06'
call KillEnterData(ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR07'
else
write(*,*) 'Block <BeginProfile>, <EndProfile> not found'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR08'
endif BF
else if (iflag == 0) then flag0
!Get the number of values
call GetData(NDEPTHS, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'NDEPTHS', &
ClientModule = 'FillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR09'
!Allocates auxiliar variables
allocate (Values (NDEPTHS))
allocate (Depth (NDEPTHS))
!Gets Depth
call GetData (Depth, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'DEPTH_PROFILE', &
ClientModule = 'FillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR10'
if (iflag /= NDEPTHS) then
write (*,*) 'Invalid Number of Depth Values. Keyword DEPTH_PROFILE'
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR11'
endif
!Gets Profile
call GetData (Values, Me%ObjEnterData, iflag, &
SearchType = ExtractType, &
keyword = 'PROFILE_VALUES', &
ClientModule = 'FillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR12'
if (iflag /= NDEPTHS) then
write (*,*) 'Invalid Number of Values. Keyword PROFILE_VALUES'
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR13'
endif
endif flag0
do k = KLB, KUB
do j = JLB, JUB
do i = ILB, IUB
if (PointsToFill3D(i, j, k) == WaterPoint) then
CellDepth = (SZZ(i, j, k) + SZZ(i, j, k - 1)) / 2 - SZZ(i, j, KUB)
Me%Matrix3D(i,j,k) = InterpolateProfile (CellDepth, NDEPTHS, Depth, Values)
endif
enddo
enddo
enddo
!Deallocates auxiliar variables
deallocate (Values)
deallocate (Depth )
call UnGetGeometry (Me%ObjGeometry, SZZ, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructSpaceProfile - ModuleFillMatrix - ERR14'
end subroutine ConstructSpaceProfile
!--------------------------------------------------------------------------
subroutine ConstructProfileTimeSerie (PointsToFill3D, ExtractType)
!Arguments-------------------------------------------------------------
integer :: ExtractType
integer, dimension(:, :, :), pointer :: PointsToFill3D
!Local----------------------------------------------------------------
real, dimension(:,:,:), pointer :: SZZ
type(T_Time) :: AuxTime
character(PathLength) :: FileName
integer :: STAT_CALL
integer :: iflag
integer :: ILB, IUB, JLB, JUB, KLB, KUB
integer :: ObjEnterData = 0
integer :: ClientNumber
integer :: FirstLine, LastLine
logical :: BlockFound
logical :: FoundSecondInstant
integer :: line, l
real, dimension(:), allocatable :: Aux
type(T_Time) :: NextTime, PreviousTime, Now
type(T_Time) :: LastInstantTime, EndTime
!Begin----------------------------------------------------------------
ILB = Me%Size3D%ILB
IUB = Me%Size3D%IUB
JLB = Me%Size3D%JLB
JUB = Me%Size3D%JUB
KLB = Me%Size3D%KLB
KUB = Me%Size3D%KUB
!Gets Center of the cells
call GetGeometryDistances(Me%ObjGeometry, SZZ = SZZ, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR10'
call GetData(FileName, &
Me%ObjEnterData , iflag, &
SearchType = ExtractType, &
keyword = 'FILENAME', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR20'
flag0: if (iflag /= 0) then
!Opens File
call ConstructEnterData(ObjEnterData, FileName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR30'
call ExtractBlockFromBuffer(ObjEnterData, ClientNumber, &
BeginTimes, EndTimes, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR40'
if (BlockFound) then
Me%ProfileTimeSerie%NumberOfInstants = LastLine - FirstLine - 1
!Allocates auxiliar variables
allocate (Me%ProfileTimeSerie%TimeInstants (Me%ProfileTimeSerie%NumberOfInstants))
l = 1
do line = FirstLine + 1, LastLine - 1
call GetData(AuxTime, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = line, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR50'
Me%ProfileTimeSerie%TimeInstants (l) = AuxTime
l = l + 1
enddo
call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR60'
else
write(*,*) 'Block <BeginTime>, <EndTime> not found'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR70'
endif
call RewindBuffer(ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR80'
call ExtractBlockFromBuffer(ObjEnterData, ClientNumber, &
BeginProfileValues, EndProfileValues, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR90'
BF: if (BlockFound) then
allocate(Aux(1:Me%ProfileTimeSerie%NumberOfInstants))
call GetData(Aux, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = FirstLine+1, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR100'
Me%ProfileTimeSerie%nValues = LastLine - FirstLine - 1
if(iflag .ne. Me%ProfileTimeSerie%NumberOfInstants)then
write(*,*) 'Number of time instants is not equal to number of values'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR110'
endif
!Allocates auxiliar variables
allocate (Me%ProfileTimeSerie%Values (1:Me%ProfileTimeSerie%nValues, &
1:Me%ProfileTimeSerie%NumberOfInstants))
l = 1
do line = FirstLine + 1, LastLine - 1
call GetData(Aux, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = line, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR120'
Me%ProfileTimeSerie%Values(l,:) = Aux(:)
l = l + 1
enddo
deallocate(Aux)
call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR130'
call RewindBuffer(ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR140'
else
write(*,*) 'Block <BeginProfileValues>, <EndProfileValues> not found'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR150'
endif BF
call ExtractBlockFromBuffer(ObjEnterData, ClientNumber, &
BeginDepth, EndDepth, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR160'
if (BlockFound) then
allocate(Aux(1:Me%ProfileTimeSerie%NumberOfInstants))
Me%ProfileTimeSerie%nDepths = LastLine - FirstLine - 1
call GetData(Aux, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = FirstLine+1, STAT = STAT_CALL)
if(iflag .ne. Me%ProfileTimeSerie%NumberOfInstants)then
write(*,*) 'Number of depths is not equal to number of values'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR170'
endif
!Allocates auxiliar variables
allocate (Me%ProfileTimeSerie%Depths (1:Me%ProfileTimeSerie%nDepths, &
1:Me%ProfileTimeSerie%NumberOfInstants))
l = 1
do line = FirstLine + 1, LastLine - 1
call GetData(Aux, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = line, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR180'
Me%ProfileTimeSerie%Depths (l,:) = Aux(:)
l = l + 1
enddo
deallocate(Aux)
call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR190'
else
write(*,*) 'Block <BeginProfileValues>, <EndProfileValues> not found'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR200'
endif
call KillEnterData(ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR210'
else if (iflag == 0) then flag0
write(*,*) 'Could not find profile time serie file.'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR220'
endif flag0
call GetComputeTimeLimits(Me%ObjTime, EndTime = EndTime, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR230'
call GetComputeCurrentTime(Me%ObjTime, Now, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR240'
if(Me%ProfileTimeSerie%NumberOfInstants > 1)then
Me%ProfileTimeSerie%PreviousInstant = 1
Me%ProfileTimeSerie%NextInstant = Me%ProfileTimeSerie%PreviousInstant
PreviousTime = Me%ProfileTimeSerie%TimeInstants(Me%ProfileTimeSerie%PreviousInstant)
call CheckCyclicMonths(PreviousTime, RefTime = Now, CyclicTimeON = Me%ProfileTimeSerie%CyclicTimeON)
i23: if (Me%ProfileTimeSerie%CyclicTimeON) then
if (Me%ProfileTimeSerie%NumberOfInstants /= 12) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR245'
else i23
if(PreviousTime .gt. Now)then
write(*,*)
write(*,*)'Could not create solution from profile time serie file'
write(*,*)'First file instant greater than current time'
write(*,*)'Matrix name: '//trim(Me%PropertyID%Name)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR250'
end if
LastInstantTime = Me%ProfileTimeSerie%TimeInstants(Me%ProfileTimeSerie%NumberOfInstants)
if (Me%ProfileTimeSerie%CyclicTimeON) call CheckCyclicMonths(LastInstantTime, RefTime = EndTime)
if(LastInstantTime .lt. EndTime)then
write(*,*)
write(*,*)'Could not create solution from profile time serie file'
write(*,*)'Last instant in file lower than simulation ending time'
write(*,*)'Matrix name: '//trim(Me%PropertyID%Name)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR260'
end if
endif i23
FoundSecondInstant = .false.
!if number of instants greater than 1 then
!find first and second instants
do while(.not. FoundSecondInstant)
Me%ProfileTimeSerie%PreviousInstant = Me%ProfileTimeSerie%NextInstant
Me%ProfileTimeSerie%NextInstant = Me%ProfileTimeSerie%NextInstant + 1
if (Me%ProfileTimeSerie%CyclicTimeON .and. Me%ProfileTimeSerie%NextInstant &
.gt. Me%ProfileTimeSerie%NumberOfInstants) then
Me%ProfileTimeSerie%NextInstant = 1
end if
NextTime = Me%ProfileTimeSerie%TimeInstants(Me%ProfileTimeSerie%NextInstant)
if (Me%ProfileTimeSerie%CyclicTimeON) call CheckCyclicMonths(NextTime, PreviousTime = PreviousTime)
if(PreviousTime .le. Now .and. NextTime .ge. Now) then
FoundSecondInstant = .true.
exit
end if
PreviousTime = NextTime
if(Me%ProfileTimeSerie%NextInstant .gt. Me%ProfileTimeSerie%NumberOfInstants &
.and. .not. Me%ProfileTimeSerie%CyclicTimeON)then
write(*,*)
write(*,*)'Could not read solution from Profile Time Serie file'
write(*,*)'Could not find second instant in file'
write(*,*)'Matrix name: '//trim(Me%PropertyID%Name)
stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR270'
end if
end do
Me%ProfileTimeSerie%PreviousTime = PreviousTime
Me%ProfileTimeSerie%NextTime = NextTime
endif
call UnGetGeometry (Me%ObjGeometry, SZZ, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR280'
allocate(Me%ProfileTimeSerie%PreviousField3D (ILB:IUB, JLB:JUB, KLB:KUB))
allocate(Me%ProfileTimeSerie%NextField3D (ILB:IUB, JLB:JUB, KLB:KUB))
Me%ProfileTimeSerie%PreviousField3D(:,:,:) = FillValueReal
Me%ProfileTimeSerie%NextField3D (:,:,:) = FillValueReal
call ProfileTimeSerieField(PointsToFill3D, Me%ProfileTimeSerie%PreviousInstant, &
Me%ProfileTimeSerie%PreviousField3D)
call ProfileTimeSerieField(PointsToFill3D, Me%ProfileTimeSerie%NextInstant, &
Me%ProfileTimeSerie%NextField3D )
if (Me%PropertyID%IsAngle) then
call InterpolateAngle3DInTime (ActualTime = Now, &
Size = Me%WorkSize3D, &
Time1 = Me%ProfileTimeSerie%PreviousTime, &
Matrix1 = Me%ProfileTimeSerie%PreviousField3D, &
Time2 = Me%ProfileTimeSerie%NextTime, &
Matrix2 = Me%ProfileTimeSerie%NextField3D, &
MatrixOut = Me%Matrix3D, &
PointsToFill3D = PointsToFill3D)
else
call InterpolateMatrix3DInTime(ActualTime = Now, &
Size = Me%WorkSize3D, &
Time1 = Me%ProfileTimeSerie%PreviousTime, &
Matrix1 = Me%ProfileTimeSerie%PreviousField3D, &
Time2 = Me%ProfileTimeSerie%NextTime, &
Matrix2 = Me%ProfileTimeSerie%NextField3D, &
MatrixOut = Me%Matrix3D, &
PointsToFill3D = PointsToFill3D)
endif
end subroutine ConstructProfileTimeSerie
!--------------------------------------------------------------------------
subroutine ConstructProfileTSDefault (PointsToFill3D, ExtractType)
!Arguments-------------------------------------------------------------
integer :: ExtractType
integer, dimension(:, :, :), pointer :: PointsToFill3D
!Local----------------------------------------------------------------
real, dimension(:,:,:), pointer :: SZZ
type(T_Time) :: AuxTime
character(PathLength) :: FileName
integer :: STAT_CALL
integer :: iflag
integer :: ILB, IUB, JLB, JUB, KLB, KUB
integer :: ObjEnterData = 0
integer :: ClientNumber
integer :: FirstLine, LastLine
logical :: BlockFound
logical :: FoundSecondInstant
integer :: line, l
real, dimension(:), allocatable :: Aux
type(T_Time) :: NextTime, PreviousTime, Now
type(T_Time) :: LastInstantTime, EndTime
!Begin----------------------------------------------------------------
ILB = Me%Size3D%ILB
IUB = Me%Size3D%IUB
JLB = Me%Size3D%JLB
JUB = Me%Size3D%JUB
KLB = Me%Size3D%KLB
KUB = Me%Size3D%KUB
!Gets Center of the cells
call GetGeometryDistances(Me%ObjGeometry, SZZ = SZZ, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructProfileTSDefault - ModuleFillMatrix - ERR10'
call GetData(FileName, &
Me%ObjEnterData , iflag, &
SearchType = ExtractType, &
keyword = 'FILENAME_DEFAULT', &
ClientModule = 'ModuleFillMatrix', &
STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) &
stop 'ConstructProfileTSDefault - ModuleFillMatrix - ERR20'
flag0: if (iflag /= 0) then
!Opens File
call ConstructEnterData(ObjEnterData, FileName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructProfileTimeSerie - ModuleFillMatrix - ERR30'
call ExtractBlockFromBuffer(ObjEnterData, ClientNumber, &
BeginTimes, EndTimes, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructProfileTSDefault - ModuleFillMatrix - ERR40'
if (BlockFound) then
Me%ProfileTimeSerie%NumberOfInstants = LastLine - FirstLine - 1
!Allocates auxiliar variables
allocate (Me%ProfileTimeSerie%TimeInstants (Me%ProfileTimeSerie%NumberOfInstants))
l = 1
do line = FirstLine + 1, LastLine - 1
call GetData(AuxTime, EnterDataID = ObjEnterData, flag = iflag, &
SearchType = FromBlock, Buffer_Line = line, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructProfileTSDefault - ModuleFillMatrix - ERR50'
Me%ProfileTimeSerie%TimeInstants (l) = AuxTime
l = l + 1
enddo
call Block_Unlock(ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructProfileTSDefault - ModuleFillMatrix - ERR60'
else
write(*,*) 'Block <BeginTime>, <EndTime> not found'
write(*,*) 'FileName = ', trim(FileName)
stop 'ConstructProfileTSDefault - ModuleFillMatrix - ERR70'
endif
call RewindBuffer(ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ConstructProfileTSDefault - ModuleFillMatrix - ERR80'