Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
7608 lines (5522 sloc) 317 KB
!------------------------------------------------------------------------------
! IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE : Mohid Model
! PROJECT : Mohid Base 2
! MODULE : Geometry
! URL : http://www.mohid.com
! AFFILIATION : IST/MARETEC, Marine Modelling Group
! DATE : May 2003
! REVISION : Frank Braunschweig - v4.0
! DESCRIPTION : Module to calculate the vertical Geometry
!
!-----------------------------------------------------------------------------
! IMPERMEABILITY : 0/1 - !Consider impermeable cell faces
! IMPER_COEF_U : real [1] !
! IMPER_COEFX_U : real [0] !
! IMPER_COEF_V : real [1] !
! IMPER_COEFX_V : real [0] !
!<begindomain>
! ID : int - !Domain ID
! TYPE : char - !Type of vertical coordinate of the domain
! !Multiple options: FIXSPACING, SIGMA,
! !CARTESIAN, FIXSEDIMENT
! LAYERS : int - !Number of layers
! EQUIDISTANT : real [0] !Equidistant layers spacing in meters
! LAYERTHICKNESS : real vector - !If not equidistant specifies layers thickness
! !starting from bottom layer (e.g. 50. 20. 10. 5.)
! TOLERANCEDEPTH : real [0.05] !Just for SIGMA,ISOPYCNIC coordinates
! TOTALTHICKNESS : real - !Total domain thickness
! !(Just for FIXSPACING, FIXSEDIMENT, SOIL_TOPLAYER)
! EMPTY_TOP_LAYERS : int [0] !Number of empty layers counting from top
! DOMAINDEPTH : real
! LAGRANGIAN : 0/1 [0] !Use lagrangian approach for distorting grometry?
! !Layers are displaced with vertical velocity
! MINEVOLVELAYERTHICKNESS : real [0.5] !Allowed distortion in percentage of initial thickness
! !(if LAGRANGIAN : 1)
! DISPLACEMENT_LIMIT : real [1000] !Maximum displacement in meters (if LAGRANGIAN : 1)
!------------------------------------------------------------------------------
!
!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 ModuleGeometry
use ModuleGlobalData
use ModuleTime
use ModuleEnterData
use ModuleGridData, only: GetGridData, UngetGridData, GetMaximumValue, &
GetGridDataFileName, WriteGridData, &
GetGridDataEvolution
use ModuleHorizontalMap, only: GetWaterPoints2D, GetExteriorBoundaryFaces, &
UnGetHorizontalMap
use ModuleHorizontalGrid, only: GetHorizontalGridSize, GetHorizontalGrid, &
GetGridCoordType, GetGridAngle, GetGridZone, &
GetLatitudeLongitude, GetGridOrigin, &
GetGridLatitudeLongitude, GetCoordTypeList, &
GetCheckDistortion, UnGetHorizontalGrid, &
GetDDecompWorkSize2D, GetDDecompParameters
#ifdef _USE_MPI
use ModuleHorizontalGrid, only: ReceiveSendLogicalMPI, GetKfloorZminMPI
#endif _USE_MPI
use ModuleFunctions, only: SetMatrixValue, SetMatrixValueAllocatable, &
Chunk_J, Chunk_K, GetPointer, ComputeAvgVerticalVelocity
use ModuleHDF5
use ModuleStopWatch, only : StartWatch, StopWatch
implicit none
private
!Subroutine----------------------------------------------------------------
!Constructor
public :: ConstructGeometry
private :: AllocateInstance
private :: ConstructGlobalVariables
private :: GetDomainsFromFile
private :: Add_Domain
private :: ComputeLayers
private :: AllocateVariables
private :: VerifyBathymetry
private :: ConstructKFloor
private :: RemoveLandBottomLayers
private :: Reallocate3DGeometryProp
private :: Reallocate3DGeometryPropR8
#ifdef _USE_SEQASSIMILATION
!Subroutine to point to memory space of geometry variables
public :: PointToGeometryState
#endif _USE_SEQASSIMILATION
!Modifier
public :: ComputeInitialGeometry !To call when the initial SurfaceElevation is known
public :: ComputeVerticalGeometry !To call in the "working cycle"
public :: UpdateKfloor
private :: ComputeSZZ
private :: ComputeFixSpacing
! private :: ComputeFixSediment
private :: ComputeInitSediment
private :: ComputeCartesian
! private :: ComputeHarmonic
private :: ComputeSigma
private :: ComputeLagrangianNew
private :: ComputeZCellCenter
private :: ComputeDistances
private :: ComputeAreas
private :: ComputeVolumes
private :: StoreVolumeZOld
private :: ComputeVolume2D
#ifdef _USE_SEQASSIMILATION
!Copy subroutines usable in sequential data assimilation to change variables' value
public :: CopyGeometryDistances
public :: CopyGeometryAreas
public :: CopyGeometryVolumes
public :: CopyGeometryWaterColumn
#endif _USE_SEQASSIMILATION
!Input / Output
public :: ReadGeometryBin
public :: WriteGeometryBin
public :: ReadGeometryHDF
public :: WriteGeometryHDF
!Selector
public :: GetGeometryDistances !SZZ, DZZ, DWZ, ZCellCenter, DUZ, DVZ, DWZ_Xgrad, DWZ_Ygrad,
public :: GetGeometryAreas
public :: GetGeometryVolumes
public :: GetGeometryKFloor
public :: GetGeometryKTop
public :: GetGeometryWaterColumn !HT
public :: GetGeometryInputFile
public :: GetGeometrySize
public :: GetGeometryMinWaterColumn
public :: GetLayer4Level
public :: UnGetGeometry
#ifdef _USE_SEQASSIMILATION
!Set subroutines usable to point variables to external variables/memory space
public :: SetGeometryDistances
public :: SetGeometryAreas
public :: SetGeometryVolumes
public :: SetGeometryWaterColumn
!Reset subroutine usable to reestablish variables to internal memory space
public :: ReSetGeometry
#endif _USE_SEQASSIMILATION
!Destructor
public :: KillGeometry
private :: DeallocateInstance
private :: DeallocateVariables
#ifdef _USE_SEQASSIMILATION
!Subroutine to point to memory space of geometry properties
public :: NullifyGeometryStatePointer
#endif _USE_SEQASSIMILATION
!Management
private :: Ready
private :: LocateObjGeometry
!Interfaces----------------------------------------------------------------
private :: UnGetGeometry2Dreal4
private :: UnGetGeometry2Dreal8
private :: UnGetGeometry3Dreal4
private :: UnGetGeometry3Dreal8
private :: UnGetGeometry2Dinteger
interface UnGetGeometry
module procedure UnGetGeometry2Dreal4
module procedure UnGetGeometry2Dreal8
module procedure UnGetGeometry3Dreal4
module procedure UnGetGeometry3Dreal8
module procedure UnGetGeometry2Dinteger
end interface UnGetGeometry
private :: ConstructGeometryV1
private :: ConstructGeometryV2
interface ConstructGeometry
module procedure ConstructGeometryV1
module procedure ConstructGeometryV2
end interface
!Parameter-----------------------------------------------------------------
!STAT
real , parameter :: FillValueDouble = -9.9e15
!Domain types
integer, parameter :: FixSpacing = 1
integer, parameter :: Sigma = 2
integer, parameter :: Isopycnic = 3
!integer, parameter :: Lagrangian = 4
integer, parameter :: Cartesian = 5
!integer, parameter :: Harmonic = 6
integer, parameter :: FixSediment = 7
integer, parameter :: SigmaTop = 8
integer, parameter :: CartesianTop = 9
!Other
integer, parameter :: INITIALGEOMETRY = 1
integer, parameter :: TRANSIENTGEOMETRY = 2
!Faces Option
integer, parameter :: MinTickness = 3
integer, parameter :: AverageTickness = 2
!Type----------------------------------------------------------------------
type T_Domain
integer :: ID = FillValueInt
integer :: DomainType = FillValueInt
logical :: IsLagrangian = .false. !Lagrangian change w/ vert veloc
integer :: NumberOfLayers = FillValueInt
integer :: EmptyTopLayers = FillValueInt
integer :: UpperLayer, LowerLayer = FillValueInt
integer :: ActiveUpperLayer = FillValueInt
real :: DomainDepth = FillValueReal
real :: TotalThickness = FillValueReal
real, dimension(:), pointer :: LayerThickness => null()
real, dimension(:), pointer :: LayerMinThickness => null()
real, dimension(:), pointer :: LayerMaxThickness => null()
real :: ToleranceDepth = FillValueReal
real :: MinInitialLayerThickness = FillValueReal
real :: MaxThicknessGrad = FillValueReal
real :: MinEvolveLayerThickness = FillValueReal
real :: MinEsp = FillValueReal
real :: BottomLayerThickness = FillValueReal
real :: GridMovementDump = FillValueReal
real :: DisplacementLimit = FillValueReal
real :: RelaxToAverageFactor = FillValueReal
integer :: InitializationMethod = FillValueInt
real :: Equidistant = FillValueReal
logical :: RomsDistortion = .false.
real :: theta_s = null_real, &
theta_b = null_real, &
Hc = null_real
logical :: SigmaZleveHybrid = .false.
real :: SigmaZleveHybrid_Hmin = null_real
real :: SigmaZleveHybrid_Hmax = null_real
type (T_Domain), pointer :: Next => null(), &
Prev => null()
end type T_Domain
type T_Distances
real, dimension(:, :, :), allocatable :: SZZ
real, dimension(:, :, :), allocatable :: DZZ
real, dimension(:, :, :), allocatable :: DWZ, DUZ, DVZ, DZI, DZE, DWZ_Xgrad, DWZ_Ygrad
real, dimension(:, :, :), allocatable :: InitialSZZ
real, dimension(:, :, :), allocatable :: ZCellCenter !Distance from the refernce level,
! center of cells, positive upwards
end type T_Distances
type T_Areas
real, dimension(:, :, :), allocatable :: AreaU, AreaV
logical :: Impermeability = .false.
real, dimension( :), allocatable :: Coef_U, CoefX_U, Coef_V, CoefX_V
end type T_Areas
type T_Volumes
real(8), dimension(:, :, :), allocatable :: VolumeZ, VolumeU, VolumeV, VolumeW, VolumeZOld
real(8), dimension(:, :), allocatable :: VolumeZ_2D
logical :: FirstVolW = .true.
end type T_Volumes
type T_KFloor
integer, dimension(:, :), allocatable :: Z, U, V, Domain
end type T_KFloor
type T_KTop
integer, dimension(:, :), allocatable :: Z
end type T_KTop
type T_WaterColumn
real, dimension(:, :), pointer :: U => null(), &
V => null(), &
Z => null() !Former HT
real :: Zmin = FillValueReal !Former Hmin
end type T_WaterColumn
#ifdef _USE_SEQASSIMILATION
type T_StatePointer
real, dimension(:, :, :), pointer :: SZZ => null(), &
DWZ => null(), &
DUZ => null(), &
DVZ => null(), &
DZZ => null(), &
ZCellCenter => null(), &
AreaU => null(), &
AreaV => null()
real, dimension(:, :), pointer :: WaterColumnU => null(), &
WaterColumnV => null(), &
WaterColumnZ => null()
real(8), dimension(:, :, :), pointer :: VolumeZ => null(), &
VolumeU => null(), &
VolumeV => null(), &
VolumeZOld => null()
end type T_StatePointer
#endif _USE_SEQASSIMILATION
type T_External
logical :: ContinuesCompute = .false.
logical :: NonHydrostatic = .false.
real :: BathymTopoFactor = 1.0
logical :: StopOnBathymetryChange = .true.
real, dimension(:,:,:), pointer :: DecayTime => null()
end type T_External
type T_Geometry
integer :: InstanceID = null_int !initialization: Jauch
integer :: FacesOption = null_int !initialization: Jauch
type (T_External) :: ExternalVar
type (T_Distances) :: Distances
type (T_Areas) :: Areas
type (T_Volumes) :: Volumes
type (T_WaterColumn) :: WaterColumn
type (T_Domain), pointer :: FirstDomain
type (T_Domain), pointer :: LastDomain
type (T_KFloor) :: KFloor
type (T_KTop) :: KTop
type (T_Time) :: ActualTime
type (T_Size3D) :: Size
type (T_Size3D) :: WorkSize
logical :: IsWindow = .false. !initialization: Jauch
logical :: LagrangianLimitsComputed = .false.
logical :: BathymNotCorrect = .false.
character(len=Pathlength) :: InputFile = null_str !initialization: Jauch
real, dimension(:,:,:), allocatable :: NearbyAvgVel_Z ! Joao Sobrinho
logical :: RemoveLandBotLayersON
#ifdef _USE_SEQASSIMILATION
!This variable is used to retain location of original memory space for variables
!changed in sequential data assimilation (some external memory is used ocasionally)
type(T_StatePointer ) :: AuxPointer
#endif _USE_SEQASSIMILATION
!Instance of other modules
integer :: ObjTopography = 0
integer :: ObjHorizontalGrid = 0
integer :: ObjHorizontalMap = 0
type (T_Geometry), pointer :: Next => null()
end type T_Geometry
!Global Module Variables
type (T_Geometry), pointer :: FirstGeometry => null()
type (T_Geometry), pointer :: Me => null()
contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine ConstructGeometryV1(GeometryID, GridDataID, HorizontalGridID, &
HorizontalMapID, ActualTime, &
NewDomain, SurfaceElevation, BathymTopoFactor, &
STAT, StopOnBathymetryChange)
!Arguments-------------------------------------------------------------
integer :: GeometryID
integer :: GridDataID
integer :: HorizontalGridID
integer :: HorizontalMapID
type (T_Time), optional :: ActualTime
character (len=*), intent(IN), optional :: NewDomain
real, dimension(:,:), pointer, optional :: SurfaceElevation
real, intent(IN), optional :: BathymTopoFactor
logical, intent(IN), optional :: StopOnBathymetryChange
integer, intent(OUT), optional :: STAT
!Local-----------------------------------------------------------------
integer :: STAT_, ready_
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mGeometry_)) then
nullify (FirstGeometry)
call RegisterModule (mGeometry_)
endif
call Ready(GeometryID, ready_)
if (ready_ .EQ. OFF_ERR_) then
!Allocates Instance
call AllocateInstance
!Associates External Instances
Me%ObjTopography = AssociateInstance (mGRIDDATA_, GridDataID )
Me%ObjHorizontalMap = AssociateInstance (mHORIZONTALMAP_, HorizontalMapID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)
!Actualize the time
if (present(ActualTime)) then
Me%ActualTime = ActualTime
endif
!Factor For MOHID Water / MOHID Land
if (present(BathymTopoFactor)) then
Me%ExternalVar%BathymTopoFactor = BathymTopoFactor
else
Me%ExternalVar%BathymTopoFactor = 1.0
endif
!stop the model on bathymetry change? - it will be verified in ConstructGlobalVariables and VerifyBathymetry
if (present(StopOnBathymetryChange)) then
Me%ExternalVar%StopOnBathymetryChange = StopOnBathymetryChange
else
Me%ExternalVar%StopOnBathymetryChange = .true.
endif
!Construct the variable common to all module
if (present(NewDomain)) then
if (present(SurfaceElevation)) then
call ConstructGlobalVariables(NewDomain = NewDomain, SurfaceElevation = SurfaceElevation)
else
call ConstructGlobalVariables(NewDomain = NewDomain)
endif
else
if (present(SurfaceElevation)) then
call ConstructGlobalVariables (SurfaceElevation = SurfaceElevation)
else
call ConstructGlobalVariables
endif
endif
!Constructs the Matrixes which contains the Ks - KFloor, etc...
if (present(SurfaceElevation)) then
call ConstructKFloor (SurfaceElevation)
else
call ConstructKFloor
endif
!Returns ID
GeometryID = Me%InstanceID
STAT_ = SUCCESS_
else
stop 'Geometry - ConstructGeometryV1 - ERR99'
end if
if (present(STAT)) STAT = STAT_
end subroutine ConstructGeometryV1
!--------------------------------------------------------------------------
subroutine ConstructGeometryV2(GeometryID, GridDataID, HorizontalGridID, &
HorizontalMapID, KMAX, &
STAT)
!Arguments-------------------------------------------------------------
integer :: GeometryID
integer :: GridDataID
integer :: HorizontalGridID
integer :: HorizontalMapID
integer :: KMAX
integer, intent(OUT), optional :: STAT
!Local-----------------------------------------------------------------
type (T_Size2D) :: WorkSize2D, Size2D
integer :: STAT_, ready_, STAT_CALL
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mGeometry_)) then
nullify (FirstGeometry)
call RegisterModule (mGeometry_)
endif
call Ready(GeometryID, ready_)
if (ready_ .EQ. OFF_ERR_) then
!Allocates Instance
call AllocateInstance
!Associates External Instances
Me%ObjTopography = AssociateInstance (mGRIDDATA_, GridDataID )
Me%ObjHorizontalMap = AssociateInstance (mHORIZONTALMAP_, HorizontalMapID )
Me%ObjHorizontalGrid = AssociateInstance (mHORIZONTALGRID_, HorizontalGridID)
!Gets horizontal size from the Bathymetry
call GetHorizontalGridSize(Me%ObjHorizontalGrid, &
Size = Size2D, &
WorkSize = WorkSize2D, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGeometryV2 - Geometry - ERR10'
Me%Size%ILB = Size2D%ILB
Me%Size%IUB = Size2D%IUB
Me%Size%JLB = Size2D%JLB
Me%Size%JUB = Size2D%JUB
Me%WorkSize%ILB = WorkSize2D%ILB
Me%WorkSize%IUB = WorkSize2D%IUB
Me%WorkSize%JLB = WorkSize2D%JLB
Me%WorkSize%JUB = WorkSize2D%JUB
call AllocateVariables(Kmax)
!Returns ID
GeometryID = Me%InstanceID
STAT_ = SUCCESS_
else
stop 'Geometry - ConstructGeometryV2 - ERR99'
end if
if (present(STAT)) STAT = STAT_
end subroutine ConstructGeometryV2
!--------------------------------------------------------------------------
subroutine AllocateInstance
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
type (T_Geometry), pointer :: NewObjGeometry
type (T_Geometry), pointer :: PreviousObjGeometry
!Allocates new instance
allocate (NewObjGeometry)
nullify (NewObjGeometry%Next)
!Insert New Instance into list and makes Current point to it
if (.not. associated(FirstGeometry)) then
FirstGeometry => NewObjGeometry
Me => NewObjGeometry
else
PreviousObjGeometry => FirstGeometry
Me => FirstGeometry%Next
do while (associated(Me))
PreviousObjGeometry => Me
Me => Me%Next
enddo
Me => NewObjGeometry
PreviousObjGeometry%Next => NewObjGeometry
endif
Me%InstanceID = RegisterNewInstance (mGEOMETRY_)
end subroutine AllocateInstance
!--------------------------------------------------------------------------
subroutine ConstructGlobalVariables(NewDomain, SurfaceElevation)
!Arguments-------------------------------------------------------------
character (len=*), intent(IN), optional :: NewDomain
real, dimension(:,:), pointer, optional :: SurfaceElevation
!External--------------------------------------------------------------
integer :: STAT_CALL
character(len = StringLength) :: Message
type (T_Size2D) :: Size2D, WorkSize2D
!----------------------------------------------------------------------
!Nullify ObjGeometryVariables
nullify (Me%FirstDomain)
!Nullify T_Volumes
! nullify (Me%Volumes%VolumeZ)
! nullify (Me%Volumes%VolumeU)
! nullify (Me%Volumes%VolumeV)
! nullify (Me%Volumes%VolumeW)
! nullify (Me%Volumes%VolumeZOld)
!Nullify T_Areas
! nullify (Me%Areas%AreaU)
! nullify (Me%Areas%AreaV)
!Nullify T_Distances
! nullify (Me%Distances%SZZ )
! nullify (Me%Distances%DZZ )
! nullify (Me%Distances%DWZ )
! nullify (Me%Distances%DUZ )
! nullify (Me%Distances%DZI )
! nullify (Me%Distances%DZE )
! nullify (Me%Distances%DVZ )
! nullify (Me%Distances%InitialSZZ )
! nullify (Me%Distances%ZCellCenter)
! nullify (Me%Distances%DWZ_Xgrad )
! nullify (Me%Distances%DWZ_Ygrad )
!Nullify T_KFloor
! nullify (Me%KFloor%Z)
! nullify (Me%KFloor%U)
! nullify (Me%KFloor%V)
! nullify (Me%KFloor%Domain)
!Gets horizontal size from the Bathymetry
call GetHorizontalGridSize(Me%ObjHorizontalGrid, &
Size = Size2D, &
WorkSize = WorkSize2D, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - Geometry - ERR02'
Me%Size%ILB = Size2D%ILB
Me%Size%IUB = Size2D%IUB
Me%Size%JLB = Size2D%JLB
Me%Size%JUB = Size2D%JUB
Me%WorkSize%ILB = WorkSize2D%ILB
Me%WorkSize%IUB = WorkSize2D%IUB
Me%WorkSize%JLB = WorkSize2D%JLB
Me%WorkSize%JUB = WorkSize2D%JUB
if (present(NewDomain)) then
Me%InputFile = NewDomain
else
!Reads the file name of the domain properties
Message ='File of the domain properties.'
call ReadFileName('DOMAIN', Me%InputFile, Message = Message, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - Geometry - ERR03'
endif
!Reads the Domain Properties
call GetDomainsFromFile
!Calculates / Verifies Layers / Computes the bounds in K
call ComputeLayers
call ConstructImpermeability
!Checks if Bathymetry is consistent with the tolerance depth - Fromer REBAIXA
!and if changes bathymetry if cartasian domain type exists
if (present(SurfaceElevation)) then
call VerifyBathymetry (SurfaceElevation)
else
call VerifyBathymetry
endif
!Allocates variables
call AllocateVariables
end subroutine ConstructGlobalVariables
!--------------------------------------------------------------------------
subroutine ConstructImpermeability
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: ObjEnterData = 0
integer :: STATUS
integer :: iflag
integer :: WorkKLB, WorkKUB
!Begin-----------------------------------------------------------------
call ConstructEnterData(ObjEnterData, Me%InputFile, STAT = STATUS)
if (STATUS /= SUCCESS_) stop "ConstructImpermeability - Geometry - ERR01"
!Searches for the FacesOption
call GetData(Me%Areas%ImPermeability, &
ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'FACES_IMPERMEABILITY', &
Default = .false., &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (Me%Areas%ImPermeability) then
allocate (Me%Areas%Coef_U (Me%Size%KLB:Me%Size%KUB))
allocate (Me%Areas%CoefX_U(Me%Size%KLB:Me%Size%KUB))
allocate (Me%Areas%Coef_V (Me%Size%KLB:Me%Size%KUB))
allocate (Me%Areas%CoefX_V(Me%Size%KLB:Me%Size%KUB))
WorkKLB = Me%WorkSize%KLB
WorkKUB = Me%WorkSize%KUB
call GetData(Me%Areas%Coef_U(WorkKLB:WorkKUB), &
ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'IMPER_COEF_U', &
Default = 1., &
ClientModule ='ModuleGeometry', &
STAT = STATUS)
call GetData(Me%Areas%CoefX_U(WorkKLB:WorkKUB), &
ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'IMPER_COEFX_U', &
Default = 0., &
ClientModule ='ModuleGeometry', &
STAT = STATUS)
call GetData(Me%Areas%Coef_V(WorkKLB:WorkKUB), &
ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'IMPER_COEF_V', &
Default = 1., &
ClientModule ='ModuleGeometry', &
STAT = STATUS)
call GetData(Me%Areas%CoefX_V(WorkKLB:WorkKUB), &
ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'IMPER_COEFX_V', &
Default = 0., &
ClientModule ='ModuleGeometry', &
STAT = STATUS)
endif
call KillEnterData(ObjEnterData, STAT = STATUS)
end subroutine ConstructImpermeability
!--------------------------------------------------------------------------
subroutine AllocateVariables(Kmax)
!Parameter-------------------------------------------------------------
integer, optional :: Kmax
!Local-----------------------------------------------------------------
integer :: STATUS
integer :: ILB, IUB, JLB, JUB, KLB, KUB
if (present(KMAX)) then
Me%WorkSize%KLB = 1
Me%WorkSize%KUB = Kmax
Me%Size%KLB = Me%WorkSize%KLB - 1
Me%Size%KUB = Me%WorkSize%KUB + 1
endif
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
KLB = Me%Size%KLB
KUB = Me%Size%KUB
!Allocates T_Volumes
allocate (Me%Volumes%VolumeZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR10'
Me%Volumes%VolumeZ = FillValueDouble
allocate (Me%Volumes%VolumeZ_2D(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR11'
Me%Volumes%VolumeZ_2D = FillValueDouble
allocate (Me%Volumes%VolumeU(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR20'
Me%Volumes%VolumeU = FillValueDouble
allocate (Me%Volumes%VolumeV(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR30'
Me%Volumes%VolumeV = FillValueDouble
allocate (Me%Volumes%VolumeZOld(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR40'
Me%Volumes%VolumeZOld = FillValueDouble
!Allocate T_Areas
allocate (Me%Areas%AreaU(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR50'
Me%Areas%AreaU = FillValueReal
allocate (Me%Areas%AreaV(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR60'
Me%Areas%AreaV = FillValueReal
!Allocate T_Distances
allocate (Me%Distances%SZZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR70'
Me%Distances%SZZ = FillValueReal
allocate (Me%Distances%DZZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR80'
Me%Distances%DZZ = FillValueReal
allocate (Me%Distances%DWZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR90'
Me%Distances%DWZ = FillValueReal
allocate (Me%Distances%DUZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR100'
Me%Distances%DUZ = FillValueReal
allocate (Me%Distances%DVZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR110'
Me%Distances%DVZ = FillValueReal
allocate (Me%Distances%InitialSZZ(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR120'
Me%Distances%SZZ = FillValueReal
allocate (Me%Distances%ZCellCenter(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR130'
Me%Distances%ZCellCenter = FillValueReal
allocate (Me%Distances%DZI(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR140'
Me%Distances%DZI = FillValueReal
allocate (Me%Distances%DZE(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR150'
Me%Distances%DZE = FillValueReal
allocate (Me%Distances%DWZ_Xgrad(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR160'
Me%Distances%DWZ_Xgrad = FillValueReal
allocate (Me%Distances%DWZ_Ygrad(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR170'
Me%Distances%DWZ_Ygrad = FillValueReal
!Allocate T_KFloor
allocate (Me%KFloor%Z(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR180'
Me%KFloor%Z = FillValueInt
allocate (Me%KFloor%U(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR190'
Me%KFloor%U = FillValueInt
allocate (Me%KFloor%V(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR200'
Me%KFloor%V = FillValueInt
allocate (Me%KFloor%Domain(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR210'
Me%KFloor%Domain = FillValueInt
allocate (Me%WaterColumn%Z(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR220'
Me%WaterColumn%Z(:,:) = FillValueReal
allocate (Me%WaterColumn%U(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR230'
Me%WaterColumn%U(:,:) = FillValueReal
allocate (Me%WaterColumn%V(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR240'
Me%WaterColumn%V(:,:) = FillValueReal
allocate (Me%KTop%Z(ILB:IUB, JLB:JUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'AllocateVariables - Geometry - ERR250'
Me%KTop%Z(:,:) = FillValueInt
allocate (Me%NearbyAvgVel_Z(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS) !Joao Sobrinho
if (STATUS /= SUCCESS_) stop "AllocateVariables - Geometry - ERR255"
Me%NearbyAvgVel_Z(:, :, :) = FillValueInt
end subroutine AllocateVariables
!--------------------------------------------------------------------------
subroutine GetDomainsFromFile
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: STATUS
integer :: FromBlock, FromFile
logical :: BlockFound
logical :: BlockLayersFound
integer :: iflag, NLayer
integer :: LBo, UBo
integer :: ClientNumber
character(len = StringLength), parameter :: block_begin = '<begindomain>'
character(len = StringLength), parameter :: block_end = '<enddomain>'
character(len = StringLength), parameter :: beginlayers = '<<beginlayers>>'
character(len = StringLength), parameter :: endlayers = '<<endlayers>>'
integer :: LagrangianOld_flag
character(len = StringLength) :: DomainType
real, dimension(:), allocatable :: AuxVector
type (T_Domain), pointer :: NewDomain
integer :: ObjEnterData = 0
integer :: i, ActualID, ID, LastLine, FirstLine, Line
LagrangianOld_flag = 0
!Get Enter data parameter
call GetExtractType(FromBlock = FromBlock, FromFile = FromFile)
call ConstructEnterData(ObjEnterData, Me%InputFile, STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR10"
!Searches for the MinWaterColumn
call GetData(Me%IsWindow, ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'WINDOW', &
Default = .false., &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR20"
!Searches for the MinWaterColumn
call GetData(Me%WaterColumn%Zmin, ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'MINIMUMDEPTH', &
Default = 0.1, &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR20"
!Searches for the FacesOption
call GetData(Me%FacesOption, ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'FACES_OPTION', &
Default = AverageTickness, &
ClientModule ='ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR30"
if (Me%FacesOption /= AverageTickness .and. Me%FacesOption /= MinTickness) then
write(*,*) "The option FACES_OPTION : ", int(Me%FacesOption), " is not valid"
write(*,*) "The only valid options are FACES_OPTION : ", AverageTickness, " advised option"
write(*,*) "The only valid options are FACES_OPTION : ", MinTickness, " advanced user option"
stop "GetDomainsFromFile - Geometry - ERR40"
endif
call GetData(Me%FacesOption, ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'FACES_OPTION', &
Default = AverageTickness, &
ClientModule ='ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR30"
call GetData(Me%RemoveLandBotLayersON, ObjEnterData, iflag, &
SearchType = FromFile, &
keyword = 'REMOVE_LAND_BOTTOM_LAYERS', &
Default = .false., &
ClientModule ='ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR40"
!The MinTickness option gives bad results espeially in the large scales.
! At the estuary scale bad results have also been identified.
!It is necessary in the future to understand the reason of this bad results.
!Rewinds Buffer
call RewindBuffer(ObjEnterData, STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR50"
!Read the defined domains
NLayer = 0
ActualID = 1
BlockFound = .true.
do while (BlockFound)
call ExtractBlockFromBuffer(ObjEnterData, ClientNumber, &
block_begin, block_end, BlockFound, &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR60"
Block: if (BlockFound) then
!Searches for the ID of the domain
call GetData(ID, ObjEnterData, iflag, &
keyword = 'ID', &
SearchType = FromBlock, &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag == 0) &
stop "GetDomainsFromFile - Geometry - ERR70"
CorretID: if (ID == ActualID) then
allocate (NewDomain)
nullify (NewDomain%Next)
nullify (NewDomain%Prev)
NewDomain%ID = ActualID
!Searches for the type of the domain
call GetData(DomainType, ObjEnterData, iflag, &
keyword = 'TYPE', &
SearchType = FromBlock, &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag == 0) &
stop "GetDomainsFromFile - Geometry - ERR80"
select case (trim(adjustl(DomainType)))
case ("FIXSPACING", "Fixspacing", "fixspacing")
NewDomain%DomainType = FixSpacing
case ("SIGMA", "Sigma", "sigma")
NewDomain%DomainType = Sigma
case ("LAGRANGIAN", "Lagrangian", "lagrangian")
!NewDomain%DomainType = Lagrangian
!for backward compatibility mark this domain as having lagrangian method
!the domain type itself will be taken from intialization method some lines below
LagrangianOld_flag = 1
NewDomain%IsLagrangian = .true.
write(*,*)
write(*,*)
write(*,*) '-------------------------------------------------------------'
write(*,*) 'WARNING: Lagrangian domains are deprecated'
write(*,*) 'new Keyword LAGRANGIAN : 1 in sigma or cartesian'
write(*,*) 'domain blocks is now used. However, for backward compability'
write(*,*) 'your domain will be processed with old keywords'
write(*,*) '-------------------------------------------------------------'
write(*,*)
write(*,*)
case ("CARTESIAN", "Cartesian", "cartesian")
NewDomain%DomainType = Cartesian
!case ("HARMONIC", "Harmonic", "harmonic")
! NewDomain%DomainType = Harmonic
case ("FIXSEDIMENT", "Fixsediment", "fixsediment")
NewDomain%DomainType = FixSediment
case ("SIGMATOP", "Sigmatop", "sigmatop")
NewDomain%DomainType = SigmaTop
case ("CARTESIANTOP", "Cartesiantop", "cartesiantop")
NewDomain%DomainType = CartesianTop
case default
stop "GetDomainsFromFile - Geometry - ERR90"
end select
!Searches for the number of Layers
call GetData(NewDomain%NumberOfLayers, ObjEnterData, iflag, &
keyword = 'LAYERS', &
SearchType = FromBlock, &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag == 0) &
stop "GetDomainsFromFile - Geometry - ERR100"
NLayer = NLayer + NewDomain%NumberOfLayers
!Allocates LayerThickness
nullify (NewDomain%LayerThickness)
LBo = NLayer - NewDomain%NumberOfLayers + 1 !This way is usefull
UBo = NLayer !in the ComputeSZZ routines
allocate (NewDomain%LayerThickness(LBo:UBo), STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR110"
nullify (NewDomain%LayerMinThickness)
nullify (NewDomain%LayerMaxThickness)
allocate (NewDomain%LayerMinThickness(LBo:UBo), STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR111"
allocate (NewDomain%LayerMaxThickness(LBo:UBo), STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR112"
! Allows the definition of equidistant layers, layer thickness = constant
call GetData(NewDomain%Equidistant, &
ObjEnterData, iflag, &
keyword = 'EQUIDISTANT', &
SearchType = FromBlock, &
ClientModule = 'ModuleGeometry', &
DEFAULT = 0. , &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR120"
cd0: if(iflag == 0) then
call ExtractBlockFromBlock(ObjEnterData, ClientNumber, &
beginlayers, endlayers, &
BlockLayersFound, &
FirstLine = FirstLine, &
LastLine = LastLine, &
STAT = STATUS)
cd1 : if (STATUS .EQ. SUCCESS_ ) then
cd2 : if (BlockLayersFound) then
if ((UBo - LBo + 1)/= (LastLine - FirstLine - 1)) &
stop "GetDomainsFromFile - Geometry - ERR130"
Line = FirstLine + 1
do i = UBo, LBo, -1
call GetData(NewDomain%LayerThickness(i), ObjEnterData, &
iflag, Buffer_Line = Line, STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag /= 1) &
stop "GetDomainsFromFile - Geometry - ERR140"
Line = Line + 1
enddo
else
allocate (AuxVector(1:UBo-LBo+1), STAT = STATUS)
!Searches for the Layer Thickness
call GetData(AuxVector, ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'LAYERTHICKNESS', &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag == 0) &
stop "GetDomainsFromFile - Geometry - ERR150"
do i = LBo, UBo
NewDomain%LayerThickness(i) = AuxVector(i-LBo+1)
enddo
deallocate (AuxVector)
endif cd2
else if (STATUS .EQ. BLOCK_END_ERR_) then cd1
stop "GetDomainsFromFile - Geometry - ERR160"
endif cd1
else cd0
do i = LBo, UBo
NewDomain%LayerThickness(i) = NewDomain%Equidistant
end do
endif cd0
!New keyword
!if geometry is lagrangian than SZZ is changed by vertical velocity
!and Lagragian approach can be used in any type of coordinates
if (LagrangianOld_flag == 0) then
call GetData(NewDomain%IsLagrangian, ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'LAGRANGIAN', &
ClientModule = 'ModuleGeometry', &
Default = .false., &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR161"
endif
!Searches for the Tolerance Depth (TYPE == SIGMA or ISOPYCNIC)
if ((NewDomain%DomainType == Sigma) .or. &
(NewDomain%DomainType == IsoPycnic)) then
call GetData(NewDomain%ToleranceDepth, ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'TOLERANCEDEPTH', &
ClientModule = 'ModuleGeometry', &
Default = 0.05, &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR170"
endif
!Searches for Total Thickness / Bottom Surface
if (NewDomain%DomainType == FixSpacing .or. &
NewDomain%DomainType == SigmaTop) then
call GetData(NewDomain%TotalThickness, ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'TOTALTHICKNESS', &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag == 0) &
stop "GetDomainsFromFile - Geometry - ERR180"
elseif(NewDomain%DomainType == FixSediment)then
call GetData(NewDomain%TotalThickness, ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'TOTALTHICKNESS', &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_)stop "GetDomainsFromFile - Geometry - ERR181"
if (iflag == 1)then
write(*,*)"TOTALTHICKNESS keyword is now obsolete when using FIXSEDIMENT."
write(*,*)"TOTALTHICKNESS keyword will be ignored."
write(*,*)"Sediment thickness is now given by a file. Check your options!"
stop "GetDomainsFromFile - Geometry - WRN180"
endif
call GetData(NewDomain%EmptyTopLayers, ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'EMPTY_TOP_LAYERS', &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_) &
stop "GetDomainsFromFile - Geometry - ERR190"
else
if (NewDomain%DomainType /= CartesianTop) then
call GetData(NewDomain%DomainDepth, ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'DOMAINDEPTH', &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag == 0) &
stop "GetDomainsFromFile - Geometry - ERR200"
endif
endif
!Searches for the coeficient which indicates how much a Lagrangian layer
!can deform
if (NewDomain%IsLagrangian) then
!The percentage of initial layer thickness that a layer can deform
!it may colapse to a size as (1 - MINEVOLVELAYERTHICKNESS) * InitialThickness
!and expand to (1 + MINEVOLVELAYERTHICKNESS)* InitialThickness
call GetData(NewDomain%MinEvolveLayerThickness, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'MINEVOLVELAYERTHICKNESS', &
ClientModule = 'ModuleGeometry', &
Default = 0.5, &
STAT = STATUS)
if (STATUS /= SUCCESS_ .or. iflag == 0) &
stop "GetDomainsFromFile - Geometry - ERR210"
if (NewDomain%MinEvolveLayerThickness < 0. .or. NewDomain%MinEvolveLayerThickness >= 1) then
write(*,*)
write(*,*)
write(*,*) 'MINEVOLVELAYERTHICKNESS keyword in Geometry Doamin'
write(*,*) 'can not be < 0 or >= 1. This keyword represents'
write(*,*) 'a percentage of initial defined thickness and the'
write(*,*) 'limit for geometry variation'
write(*,*) 'Please verify values'
write(*,*)
stop "GetDomainsFromFile - Geometry - ERR212"
endif
! call GetData(NewDomain%GridMovementDump, &
! ObjEnterData, iflag, &
! SearchType = FromBlock, &
! keyword = 'GRIDMOVEMENTDUMP', &
! ClientModule = 'ModuleGeometry', &
! Default = 0.0, &
! STAT = STATUS)
! if (STATUS /= SUCCESS_) &
! stop "GetDomainsFromFile - Geometry - ERR220"
!maximum displacement in meters. If too high will not limit more
!than MINEVOLVELAYERTHICKNESS factor
call GetData(NewDomain%DisplacementLimit, &
ObjEnterData, iflag, &
keyword = 'DISPLACEMENT_LIMIT', &
Default = 1000., &
SearchType = FromBlock, &
ClientModule = 'ModuleGeometry', &
STAT = STATUS)
if (STATUS /= SUCCESS_) &
stop "GetDomainsFromFile - Geometry - ERR230"
!Joao Sobrinho
call GetData(NewDomain%RelaxToAverageFactor, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'RELAXTOAVERAGEFACTOR', &
ClientModule = 'ModuleGeometry', &
Default = 0.7, &
STAT = STATUS)
if (STATUS /= SUCCESS_) &
stop "GetDomainsFromFile - Geometry - ERR235"
if (LagrangianOld_flag == 1) then
call GetData(DomainType, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'INITIALIZATION_METHOD ', &
ClientModule = 'ModuleGeometry', &
Default = "SIGMA", &
STAT = STATUS)
if (STATUS /= SUCCESS_) &
stop "GetDomainsFromFile - Geometry - ERR240"
!Lagrangian is no longer a domain but a process.
!As so the domaintype here will be used as the domain type
!and the previous keywords are used to compute lagrangian method
!in any of the two types
if (DomainType .EQ. "SIGMA" ) then
!NewDomain%InitializationMethod = Sigma
NewDomain%DomainType = Sigma
elseif (DomainType .EQ. "CARTESIAN" ) then
!NewDomain%InitializationMethod = Cartesian
NewDomain%DomainType = Cartesian
else
write (*,*) "Initialization Method invalid"
stop "GetDomainsFromFile - Geometry - ERR250"
endif
endif
endif
!Searches for the coeficient which indicates the minimal thickness in %
!of the bottom cells if MinInitialLayerThickness = 1 - classic cartesian coordinates
!if < 1 - cartesian with shave cells
!In cartesian coordinates the reference layer thickness is compute for the
!maximum depth
! if (NewDomain%DomainType == Cartesian .or. &
! NewDomain%DomainType == Harmonic) then
if (NewDomain%DomainType == Cartesian) then
call GetData(NewDomain%MinInitialLayerThickness, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'MININITIALLAYERTHICKNESS', &
ClientModule = 'ModuleGeometry', &
Default = 0.2, &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR270"
call GetData(NewDomain%MaxThicknessGrad, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'MAX_THICKNESS_GRAD', &
ClientModule = 'ModuleGeometry', &
Default = 2./NewDomain%MinInitialLayerThickness, &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR275"
if (NewDomain%MaxThicknessGrad < 1.) stop "GetDomainsFromFile - Geometry - ERR277"
endif
if (NewDomain%DomainType == CartesianTop) then
NewDomain%MinInitialLayerThickness = 0.0
endif
!Seraches for the minimum thickness of colapsing cells of the Harmonic domain
! if (NewDomain%DomainType == Harmonic) then
!
! call GetData(NewDomain%MinEsp, &
! ObjEnterData, iflag, &
! SearchType = FromBlock, &
! keyword = 'MIN_TOP_THICKNESS', &
! ClientModule = 'ModuleGeometry', &
! Default = Me%WaterColumn%ZMin, &
! STAT = STATUS)
! if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR280"
!
! endif
!Seraches for the minimum thickness of bottom layer
if (NewDomain%DomainType == CartesianTop) then
call GetData(NewDomain%BottomLayerThickness, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'MIN_BOTTOM_THICKNESS', &
ClientModule = 'ModuleGeometry', &
Default = Me%WaterColumn%ZMin, &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR281"
endif
if (NewDomain%DomainType == Sigma) then
call GetData(NewDomain%RomsDistortion, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'ROMS_DISTORTION', &
ClientModule = 'ModuleGeometry', &
Default = .false., &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR300"
if (NewDomain%RomsDistortion) then
!theta_s = terrain following coodinates bottom control parameter = 5
call GetData(NewDomain%theta_s, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'THETA_S', &
ClientModule = 'ModuleGeometry', &
Default = 5., &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR310"
!theta_b = terrain following coordinates surface control parameter = 0.4
call GetData(NewDomain%theta_b, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'THETA_B', &
ClientModule = 'ModuleGeometry', &
Default = 0.4, &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR320"
!hc = s coordinate parameter critical depth = 25
call GetData(NewDomain%Hc, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'HC', &
ClientModule = 'ModuleGeometry', &
Default = 25., &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR330"
endif
endif
call SigmaZleveHybridOptions(NewDomain,ObjEnterData)
!Inserts new domain into the domain list
!call AddDomainToGeometry(NewDomain)
call Add_Domain(NewDomain)
!Next ID to search for
ActualID = ActualID + 1
!Rewinds Buffer
call Block_Unlock(ObjEnterData, ClientNumber)
call RewindBuffer(ObjEnterData, STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR390"
endif CorretID
endif Block
enddo
call Block_Unlock(ObjEnterData, ClientNumber)
call KillEnterData(ObjEnterData, STAT = STATUS)
if (STATUS /= SUCCESS_) stop "GetDomainsFromFile - Geometry - ERR400"
end subroutine GetDomainsFromFile
!--------------------------------------------------------------------------
subroutine SigmaZleveHybridOptions(NewDomain, ObjEnterData)
!Parameter-------------------------------------------------------------
type (T_Domain), pointer :: NewDomain
integer :: ObjEnterData
!Local-----------------------------------------------------------------
integer :: iflag, STAT_CALL
!Begin-----------------------------------------------------------------
call GetData(NewDomain%SigmaZleveHybrid, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'SIGMA_Z_SURFACE', &
ClientModule = 'ModuleGeometry', &
Default = .false., &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop "SigmaZleveHybridOptions - Geometry - ERR10"
i1: if (NewDomain%SigmaZleveHybrid) then
call GetData(NewDomain%SigmaZleveHybrid_Hmin, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'SIGMA_Z_SURFACE_HMIN', &
ClientModule = 'ModuleGeometry', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop "SigmaZleveHybridOptions - Geometry - ERR20"
if (iflag == 0) then
write(*,*) 'if option SIGMA_Z_SURFACE is ON need to define'
write(*,*) 'SIGMA_Z_SURFACE_HMIN'
stop
endif
call GetData(NewDomain%SigmaZleveHybrid_Hmax, &
ObjEnterData, iflag, &
SearchType = FromBlock, &
keyword = 'SIGMA_Z_SURFACE_HMAX', &
ClientModule = 'ModuleGeometry', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop "SigmaZleveHybridOptions - Geometry - ERR30"
if (iflag == 0) then
write(*,*) 'if option SIGMA_Z_SURFACE is ON need to define'
write(*,*) 'SIGMA_Z_SURFACE_HMAX'
stop
endif
endif i1
end subroutine SigmaZleveHybridOptions
!--------------------------------------------------------------------------
subroutine ComputeLayers
!Parameter-------------------------------------------------------------
!Local-----------------------------------------------------------------
type (T_Domain), pointer :: CurrentDomain
real :: RelativSum, Error, &
BottomDepth, TopDepth
real(8) :: Sum, DomainDif
integer :: iLayer, LayersBelow
!Begin-----------------------------------------------------------------
!Computes limits for K
Me%WorkSize%KUB = 0
CurrentDomain => Me%FirstDomain
do while (associated(CurrentDomain))
Me%WorkSize%KUB = Me%WorkSize%KUB + CurrentDomain%NumberOfLayers
CurrentDomain => CurrentDomain%Next
enddo
Me%WorkSize%KLB = 1
Me%Size%KLB = 0
Me%Size%KUB = Me%WorkSize%KUB + 1
LayersBelow = 0
CurrentDomain => Me%FirstDomain
do while (associated(CurrentDomain))
!Computes UpperLayer and LowerLayer
CurrentDomain%LowerLayer = LayersBelow + 1
CurrentDomain%UpperLayer = LayersBelow + CurrentDomain%NumberOfLayers
if (CurrentDomain%ID == Me%LastDomain%ID) CurrentDomain%ActiveUpperLayer = CurrentDomain%UpperLayer
LayersBelow = LayersBelow + CurrentDomain%NumberOfLayers
!Verifies and corrects relativ layer thickness
!This is just done in the case of Sigma, Fixspacing, FixSediment
!initialization
if (CurrentDomain%DomainType == Sigma .or. &
CurrentDomain%DomainType == Fixspacing .or. &
CurrentDomain%DomainType == FixSediment .or. &
CurrentDomain%DomainType == SigmaTop .or. &
CurrentDomain%DomainType == IsoPycnic) then
RelativSum = 0.
do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
RelativSum = RelativSum + CurrentDomain%LayerThickness(iLayer)
enddo
Error = abs(1.-RelativSum)
if ( Error > 0.01) then
write(*,*)'Layers incorrectly defined - Domain number ',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR10.'
!PCL
else if (Error<=0.01 .and. Error > 1.e-6) then
do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
CurrentDomain%LayerThickness(iLayer) = &
CurrentDomain%LayerThickness(iLayer) - &
(1. - RelativSum) / &
real(CurrentDomain%NumberOfLayers)
enddo
endif
else if (CurrentDomain%DomainType == Cartesian) then
Sum = 0.
do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
Sum = Sum + dble(CurrentDomain%LayerThickness(iLayer))
enddo
if (CurrentDomain%ID == Me%LastDomain%ID) then
TopDepth = 0.
else
TopDepth = CurrentDomain%DomainDepth
endif
if (CurrentDomain%ID > 1) then
BottomDepth = CurrentDomain%Prev%DomainDepth
else
call GetMaximumValue(Me%ObjTopography, BottomDepth)
endif
DomainDif = BottomDepth - TopDepth
!Error = abs(DomainDif - Sum)
if ((DomainDif - Sum)>1e-7) then
if (.not. Me%IsWindow) then
write(*,*)'Layers incorrectly defined - Domain number ',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR20.'
endif
endif
endif
!Verifies if there is a fixspacing domain, other then the lower domain
if (CurrentDomain%ID /= 1 .and. CurrentDomain%DomainType == FixSpacing) then
write(*,*)'Just the first domain can be of the type FixSpacing'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR30.'
endif
!Verifies if there is a fixsediment domain, other then the lower domain
if (CurrentDomain%ID /= 1 .and. CurrentDomain%DomainType == FixSediment) then
write(*,*)'Just the first domain can be of the type FixSediment and it cant coexist with other domains'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR40.'
endif
!Verifies if a cartesian domain is above a fixspacing domain
if (CurrentDomain%ID > 1 .and. CurrentDomain%DomainType == Cartesian .and. &
(Me%FirstDomain%DomainType == FixSpacing .or. &
Me%FirstDomain%DomainType == FixSediment)) then
write(*,*)'A Cartesian domain type cant overlay a Fixspacing domain'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR50.'
endif
! !Verifies if a harmonic domain is above a fixspacing domain
! if (CurrentDomain%ID > 1 .and. CurrentDomain%DomainType == Harmonic .and. &
! (Me%FirstDomain%DomainType == FixSpacing .or. &
! Me%FirstDomain%DomainType == FixSediment)) then
! write(*,*)'A Harmonic domain type cant overlay a Fixspacing or FixSediment domain'
! write(*,*)'Verify domain number :',CurrentDomain%ID
! stop 'ComputeLayers - Geometry - ERR60.'
! endif
!Verifies if a Lagrangian domain is above a fixspacing domain
if (CurrentDomain%ID > 1 .and. CurrentDomain%IsLagrangian .and. &
(Me%FirstDomain%DomainType == FixSpacing .or. &
Me%FirstDomain%DomainType == FixSediment)) then
write(*,*)'A Largragian domain type cant overlay a Fixspacing or FixSediment domain'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR70.'
endif
!Verifies if SigmaTop domain is topmost domain
if (CurrentDomain%DomainType == SigmaTop .and. &
CurrentDomain%ID /= Me%LastDomain%ID) then
write(*,*)'A SigmaTop domain type can only be the topmost domain'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR10.'
endif
!Verifies if there is a sigma domain below a SigmaTop domain (just one layer?)
if (CurrentDomain%DomainType == SigmaTop .and. CurrentDomain%ID == 1) then
write(*,*)'A SigmaTop domain needs a Sigma domain below'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR10.'
endif
!Verifies if there is a sigma domain below a SigmaTop domain (sigma below)
if (CurrentDomain%DomainType == SigmaTop) then
if (CurrentDomain%Prev%DomainType /= SIGMA) then
write(*,*)'A SigmaTop domain needs a Sigma domain below'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR11.'
endif
endif
!Verifies if CartesianTop domain is topmost domain
if (CurrentDomain%DomainType == CartesianTop) then
if (CurrentDomain%ID /= Me%LastDomain%ID .or. CurrentDomain%ID /= 1) then
write(*,*)'A CartesianTop domain type can only be the topmost domain'
write(*,*)'A CartesianTop domain type can only exist alone'
write(*,*)'Verify domain number :',CurrentDomain%ID
stop 'ComputeLayers - Geometry - ERR11a.'
endif
endif
!Verifies if Lagrangian method is used with SIGMA or Cartesian
if (CurrentDomain%IsLagrangian) then
if (CurrentDomain%DomainType /= Cartesian .and. CurrentDomain%DomainType /= Sigma) then
write(*,*)'Lagrangian method can be used with SIGMA'
write(*,*)'or CARTESIAN domains. Verify if other is used'
stop 'ComputeLayers - Geometry - ERR20.'
endif
endif
! !Verifies if the Minimal Layer thickness of a lagrangian domain is smalles then
! !50% of the initial layer thickness
! if (CurrentDomain%DomainType == Lagrangian) then
! Error = 1.e10
! do iLayer = CurrentDomain%LowerLayer, CurrentDomain%UpperLayer
! if (Error > CurrentDomain%LayerThickness(iLayer)) then
! Error = CurrentDomain%LayerThickness(iLayer)
! endif
! enddo
! if (CurrentDomain%MinEvolveLayerThickness > 0.50) then
! write(*,*)'The MinimalThickness of the layers should not be greater then 50% of'
! write(*,*)'of the initial layer thickness.'
! write(*,*)'Verify domain number :',CurrentDomain%ID
! write(*,*)'Layer :',iLayer
! stop 'ComputeLayers - Geometry - ERR80.'
! else if (CurrentDomain%DisplacementLimit > 0.499 * Error) then
! write(*,*)'The DisplacementLimit of the layers is higher than half the minimum LayerThickness'
! stop 'ComputeLayers - Geometry - ERR90.'
! endif
! endif
!Verifies if the inital Layer thickness of the cartesian coordinates is not small then
!then 1%
! if (CurrentDomain%DomainType == Cartesian .or. &
! CurrentDomain%DomainType == Harmonic ) then
if (CurrentDomain%DomainType == Cartesian) then
if (CurrentDomain%MinInitialLayerThickness < 0.01) then
write(*,*)'The MinimalThickness of the layers should not be smaller then 1% of'
write(*,*)'of the reference layer thickness'
write(*,*)'Verify domain number :',CurrentDomain%ID
write(*,*)'Layer :',iLayer
stop 'ComputeLayers - Geometry - ERR100.'
endif
endif
CurrentDomain => CurrentDomain%Next
enddo
end subroutine ComputeLayers
!--------------------------------------------------------------------------
!In the case that the domain is Sigma, Isopycnic the Bathymetry is tested in
!order to avoid layers with zero height.
!In the case of a cartesian domain, the Bathymetry is rearranged, in order to avoid layers to thin
subroutine VerifyBathymetry (SurfaceElevation)
!Arguments-------------------------------------------------------------
real, dimension(:, :), pointer, optional :: SurfaceElevation
!Local-----------------------------------------------------------------
real, dimension(:, :), pointer :: Bathymetry, NewBathymetry, XX_IE, YY_IE
real, dimension(:), pointer :: XX, YY
real :: GRID_ANGLE, Latitude, Longitude, XORIG, YORIG
integer :: ICOORD_TIP, Zone
integer :: GEOG, UTM, MIL_PORT, SIMPLE_GEOG, NLRD, GRID_COORD
integer :: ILB, IUB, JLB, JUB
integer :: i, j, iLayer, STAT_CALL
real :: DomainDepth, Tolerance
real :: TopDepth, DomainThickness
real :: LayerTopDepth, LayerMinBottomDepth
real :: LayerTop, LayerBottom
real :: DistToBottom, DistToTop
real :: MinimalThickness, AllmostZero_
real :: BottomDepth
character(len=StringLength) :: BathymetryFile
character(len=StringLength) :: Comment1, Comment2
logical :: WriteNewBathymetry = .false., ConvertsWaterInLand
!logical :: WriteNewBathymetry = .false., Distortion, ConvertsWaterInLand
type (T_Domain), pointer :: CurrentDomain
type (T_Size2D) :: Size2D
integer :: LengthWithoutExt
real :: T1, Taux, Tmax, d1, d2
character(len=StringLength) :: FileVersion
integer :: FileVersionInt
logical :: WriteBathymAux, MasterOrSlave, Master
!Begin-----------------------------------------------------------------
ILB = Me%WorkSize%ILB; Size2D%ILB = Me%WorkSize%ILB
IUB = Me%WorkSize%IUB; Size2D%IUB = Me%WorkSize%IUB
JLB = Me%WorkSize%JLB; Size2D%JLB = Me%WorkSize%JLB
JUB = Me%WorkSize%JUB; Size2D%JUB = Me%WorkSize%JUB
call GetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR10'
nullify (NewBathymetry)
allocate(NewBathymetry(Me%Size%ILB:Me%Size%IUB, Me%Size%JLB:Me%Size%JUB))
!Copy bathymetry array
!NewBathymetry(:,:) = Bathymetry(:,:)
call SetMatrixValue(NewBathymetry, Size2D, Bathymetry)
CurrentDomain => Me%FirstDomain
do1: do while (associated(CurrentDomain))
cd1: if ((CurrentDomain%DomainType == Sigma) .or. &
(CurrentDomain%DomainType == IsoPycnic)) then
if (CurrentDomain%ID == Me%LastDomain%ID) exit do1 !Dont correct the upper domain
DomainDepth = CurrentDomain%DomainDepth
Tolerance = CurrentDomain%ToleranceDepth
do j = JLB, JUB
do i = ILB, IUB
if ((Bathymetry(i, j) > DomainDepth) .and. &
(Bathymetry(i, j) < DomainDepth + Tolerance)) then
write(*,*) NewBathymetry(i, j), DomainDepth, Tolerance
NewBathymetry(i, j) = DomainDepth
WriteNewBathymetry = .true.
endif
enddo
enddo
!Correction in the case of cartesian domain
! else if (CurrentDomain%DomainType == Cartesian .or. &
! CurrentDomain%DomainType == Harmonic) then
else if (CurrentDomain%DomainType == Cartesian) then
!Gets the upper and the lower depth of the domain and calculates its
!thickness
if (CurrentDomain%ID == Me%LastDomain%ID) then
TopDepth = 0.
else
TopDepth = CurrentDomain%DomainDepth
endif
if (CurrentDomain%ID > 1) then
BottomDepth = CurrentDomain%Prev%DomainDepth
else
call GetMaximumValue(Me%ObjTopography, BottomDepth)
endif
DomainThickness = BottomDepth - TopDepth
!Verifies if the depth of the Bathymetry is between the upper and the lower
!part of the layer, if so corrects the Bathymetry
!write(*,*)'Checking if Bathymetry is between LayerTopDepth and LayerMinBottomDepth'
LayerTopDepth = TopDepth
doLayer: do iLayer = CurrentDomain%UpperLayer, CurrentDomain%LowerLayer, -1
!MinimalThickness is now defined as a percentage of the initial thickness of each layer
MinimalThickness = CurrentDomain%MinInitialLayerThickness
!The LayerThickness is now given in meters, so dont multiply by DomainTickness
!Frank - Jan 2001
!LayerMinBottomDepth = LayerTopDepth + CurrentDomain%LayerThickness(iLayer) * &
! DomainThickness * MinimalThickness
LayerMinBottomDepth = LayerTopDepth + CurrentDomain%LayerThickness(iLayer) * &
MinimalThickness
doj: do j = JLB, JUB
doi: do i = ILB, IUB
! Avoid round-off errors - Hernani Sept 2004
! if ((Bathymetry(i, j) > LayerTopDepth) .and. &
! (Bathymetry(i, j) < LayerMinBottomDepth)) then
! GRiflet, NVerelst - Usar como erro de truncatura na escrita do ficheiro de batimetria
! metade da menor divisao: i.e. 5e-4 m
ConvertsWaterInLand = .false.
d1 = LayerTopDepth
d2 = LayerTopDepth + CurrentDomain%LayerThickness(iLayer)
AllmostZero_ = 5e-4
if (Bathymetry(i, j) >= d2) cycle
if (Bathymetry(i, j) <= d1) cycle
if ((Bathymetry(i, j) - LayerTopDepth > AllmostZero_) .and. &
(Bathymetry(i, j) - LayerMinBottomDepth < -AllmostZero_)) then
ConvertsWaterInLand = .true.
endif
if (.not. ConvertsWaterInLand) then
Tmax = FillValueReal
T1 = Bathymetry(i, j) - LayerTopDepth
if (T1 <=AllmostZero_) cycle
Taux = Bathymetry(i-1,j) - LayerTopDepth
if (Taux > CurrentDomain%LayerThickness(iLayer)) then
Taux = CurrentDomain%LayerThickness(iLayer)
endif
if ( Taux > Tmax) Tmax = Taux
Taux = Bathymetry(i+1,j) - LayerTopDepth
if (Taux > CurrentDomain%LayerThickness(iLayer)) then
Taux = CurrentDomain%LayerThickness(iLayer)
endif
if ( Taux > Tmax) Tmax = Taux
Taux = Bathymetry(i,j+1) - LayerTopDepth
if (Taux > CurrentDomain%LayerThickness(iLayer)) then
Taux = CurrentDomain%LayerThickness(iLayer)
endif
if ( Taux > Tmax) Tmax = Taux
Taux = Bathymetry(i,j-1) - LayerTopDepth
if (Taux > CurrentDomain%LayerThickness(iLayer)) then
Taux = CurrentDomain%LayerThickness(iLayer)
endif
if ( Taux > Tmax) Tmax = Taux
if (Tmax/T1 > CurrentDomain%MaxThicknessGrad) then
ConvertsWaterInLand = .true.
endif
endif
if (ConvertsWaterInLand) then
NewBathymetry(i, j) = LayerTopDepth
write(*,*)'i = ', i
write(*,*)'j = ', j
write(*,*)'Bathymetry = ', Bathymetry(i, j)
write(*,*)'New Bathymetry = ',NewBathymetry(i, j)
WriteNewBathymetry = .true.
endif
enddo doi
enddo doj
!
!The Layer thickness is now given in meters, so dont multiply by DomainThickness
!Frank Jan 2001
LayerTopDepth = LayerTopDepth + CurrentDomain%LayerThickness(iLayer)
enddo doLayer
elseif (CurrentDomain%DomainType == CartesianTop) then
do j = JLB, JUB
do i = ILB, IUB
if (Bathymetry(i, j) > -55.) then
TopDepth = SurfaceElevation(i, j)
BottomDepth = Bathymetry(i, j)
iLayer = CurrentDomain%UpperLayer
LayerTop = TopDepth
LayerBottom = LayerTop - CurrentDomain%LayerThickness(iLayer)
do while (iLayer >= CurrentDomain%LowerLayer)
AllmostZero_ = AllmostZeroFraction * CurrentDomain%LayerThickness(iLayer)
! if (LayerBottom - AllmostZero_ <= BottomDepth .and. LayerTop + AllmostZero_ >= BottomDepth) then
DistToBottom = BottomDepth - LayerBottom
DistToTop = LayerTop - BottomDepth
if ((BottomDepth > LayerBottom) .and. (BottomDepth < LayerTop) .and. &
(abs(DistToTop) > AllmostZero_) .and. (abs(DistToBottom) > AllmostZero_)) then
if ((iLayer .eq. CurrentDomain%UpperLayer) .or. &
(DistToBottom .le. DistToTop)) then
NewBathymetry(i, j) = LayerBottom
WriteNewBathymetry = .true.
write(*,*)'i = ', i
write(*,*)'j = ', j
write(*,*)'Bathymetry = ',Bathymetry(i, j)
write(*,*)'New Bathymetry = ',NewBathymetry(i, j)
else
NewBathymetry(i, j) = LayerTop
WriteNewBathymetry = .true.
write(*,*)'i = ', i
write(*,*)'j = ', j
write(*,*)'Bathymetry = ',Bathymetry(i, j)
write(*,*)'New Bathymetry = ',NewBathymetry(i, j)
endif
exit
else
iLayer = iLayer - 1
if (iLayer > 0) then
LayerTop = LayerBottom
LayerBottom = LayerTop - CurrentDomain%LayerThickness(iLayer)
else
exit
endif
endif
enddo
endif
enddo
enddo
endif cd1
CurrentDomain => CurrentDomain%Next
enddo do1
call GetDDecompParameters(HorizontalGridID = Me%ObjHorizontalGrid, &
MasterOrSlave = MasterOrSlave, &
Master = Master, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR20'
if (MasterOrSlave) then
WriteBathymAux = WriteNewBathymetry
#ifdef _USE_MPI
!Check if at least one of the sub-domains need to change the bathymetry
call ReceiveSendLogicalMPI(HorizontalGridID = Me%ObjHorizontalGrid, &
LogicalIn = WriteBathymAux, &
LogicalOut = WriteNewBathymetry, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR30'
#endif _USE_MPI
endif
if (WriteNewBathymetry) then
!Gets XX and YY
call GetHorizontalGrid(Me%ObjHorizontalGrid, XX = XX, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR40'
call GetHorizontalGrid(Me%ObjHorizontalGrid, YY = YY, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR50'
call GetCoordTypeList (GEOG = GEOG, UTM = UTM, MIL_PORT = MIL_PORT, &
SIMPLE_GEOG = SIMPLE_GEOG, GRID_COORD = GRID_COORD, &
NLRD = NLRD)
!Gets the type of Coordinates
call GetGridCoordType(Me%ObjHorizontalGrid, ICOORD_TIP, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR60'
if (ICOORD_TIP == SIMPLE_GEOG)then
call GetGridLatitudeLongitude(Me%ObjHorizontalGrid, &
GridLatitudeConn = YY_IE, &
GridLongitudeConn = XX_IE, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR70'
elseif(ICOORD_TIP == UTM .or. ICOORD_TIP == MIL_PORT .or. &
ICOORD_TIP == GRID_COORD .or. ICOORD_TIP == NLRD)then
!Gets XX_IE and YY_IE
call GetHorizontalGrid(Me%ObjHorizontalGrid, XX_IE = XX_IE, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR80'
call GetHorizontalGrid(Me%ObjHorizontalGrid, YY_IE = YY_IE, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR90'
else
write(*,*)'GEOG coordinate type cannot be used in digital terrain generation'
stop 'VerifyBathymetry - Geometry - ERR100'
end if
!Gets the Grid angle
call GetGridAngle (Me%ObjHorizontalGrid, GRID_ANGLE, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR110'
!Gets the Latitude / Longitude
call GetLatitudeLongitude(Me%ObjHorizontalGrid, Latitude, Longitude, STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR120'
!Gets Zone of the Bathymetry
call GetGridZone (Me%ObjHorizontalGrid, Zone, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR130'
!Gets Origin of the Bathymetry
call GetGridOrigin (Me%ObjHorizontalGrid, Xorig, Yorig, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR140'
!Writes new Bathymetry
call GetGridDataFileName(Me%ObjTopography, FileName = BathymetryFile, STAT= STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR150'
BathymetryFile = adjustl(BathymetryFile)
Comment1 = "Automatic Generated Grid Data File"
Comment2 = "Based On "//trim(BathymetryFile)
LengthWithoutExt= len_trim(BathymetryFile) - 4
!BathymetryFile = BathymetryFile(1:LengthWithoutExt)//"_"//".new"
!check if already has version in name (2 carachters and can go up to 99 versions)
if (BathymetryFile(LengthWithoutExt-3:LengthWithoutExt-2) == "_v") then
!read existing file version
FileVersion = trim(adjustl(BathymetryFile(LengthWithoutExt-1:LengthWithoutExt)))
read(FileVersion, *, IOSTAT = STAT_CALL) FileVersionInt
if (STAT_CALL /= 0) then
!user used a letter after "_v". it may happen in a original filename
FileVersionInt = 0
endif
!update file version
FileVersionInt = FileVersionInt + 1
write(FileVersion, "(i2)") FileVersionInt
if (FileVersionInt < 10) then
FileVersion = "0"//FileVersion(2:2)
endif
else
FileVersionInt = 1
endif
if (FileVersionInt == 1) then
BathymetryFile = BathymetryFile(1:LengthWithoutExt)//"_v"//"01.dat"
else
BathymetryFile = BathymetryFile(1:LengthWithoutExt-2)//FileVersion//".dat"
endif
call WriteGridData (FileName = BathymetryFile, &
COMENT1 = Comment1, &
COMENT2 = Comment2, &
HorizontalGridID = Me%ObjHorizontalGrid, &
FillValue = -99., &
Overwrite = ON, &
GridData2D_Real = NewBathymetry, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR160'
!call GetCheckDistortion (Me%ObjHorizontalGrid, Distortion, STAT = STAT_CALL)
!if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR140'
!if (.not. Distortion) then
! call WriteGridData (FileName = BathymetryFile, &
! XX = XX, &
! YY = YY, &
! COMENT1 = Comment1, &
! COMENT2 = Comment2, &
! WorkSize = Size2D, &
! CoordType = ICOORD_TIP, &
! Xorig = Xorig, &
! Yorig = Yorig, &
! Zone = Zone, &
! GRID_ANGLE = GRID_ANGLE, &
! Latitude = Latitude, &
! Longitude = Longitude, &
! GridData2D_Real= NewBathymetry, &
! FillValue = -99., &
! Overwrite = ON, &
! STAT = STAT_CALL)
! if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR150'
!else
! call WriteGridData (FileName = BathymetryFile, &
! ConnectionX = XX_IE, &
! ConnectionY = YY_IE, &
! COMENT1 = "***", &
! COMENT2 = "***", &
! WorkSize = Size2D, &
! CoordType = ICOORD_TIP, &
! Xorig = Xorig, &
! Yorig = Yorig, &
! Zone = Zone, &
! GRID_ANGLE = GRID_ANGLE, &
! Latitude = Latitude, &
! Longitude = Longitude, &
! GridData2D_Real= NewBathymetry, &
! FillValue = -99., &
! Overwrite = ON, &
! STAT = STAT_CALL)
! if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR160'
!endif
deallocate(NewBathymetry)
nullify (NewBathymetry)
!Displays Message to inform
!Check if will stop or not. By default always stop and ask user to run again
!But for run on demand this needs to be possible
call SetError(WARNING_, INTERNAL_, 'Bathymetry changed due to geometry', Screen = .false.)
write(*,*)'WARNING'
write(*,*)'A new Bathymetry has been created, which consists with the geometry'
write(*,*)'New Bathymetry file : ', trim(BathymetryFile)
write(*,*)''
if (Me%ExternalVar%StopOnBathymetryChange) then
if (.not.MasterOrSlave .or. (MasterOrSlave .and. Master)) then
write(*,*)'Modify the file Nomfich.dat and Re-run the model'
stop 'VerifyBathymetry - Geometry - ERR170'
endif
endif
endif
!Disposes pointer to the Bathymetry
call UngetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'VerifyBathymetry - Geometry - ERR180'
!----------------------------------------------------------------------
end subroutine VerifyBathymetry
!--------------------------------------------------------------------------
subroutine ConstructKFloor (SurfaceElevation)
!Arguments-------------------------------------------------------------
real, dimension(:, :), pointer, optional :: SurfaceElevation
!Local-----------------------------------------------------------------
real, dimension(:, :), pointer :: Bathymetry
integer, dimension(:, :), pointer :: WaterPoints2D
integer :: ILB, IUB, JLB, JUB, KUB
integer :: iLayer, i, j, LayersBelow
real(8) :: BottomDepth, TopDepth, LayerTopDepth
real(8) :: LayerBottomDepthMin, LayerBottomDepthMax
real(8) :: MinimalThickness, DomainThickness
real(8) :: LayerThicknessMin, LayerThicknessMax
logical :: FoundKFloor
type (T_Domain), pointer :: CurrentDomain
integer, dimension(:,:), pointer :: ExteriorFacesU, ExteriorFacesV
integer :: STAT_CALL
logical :: NeedToStop
real(8) :: MaxDomainThickness
real(8) :: TotalLayerThickness, AuxDepth
real :: Aux4, AllmostZero_
real(8) :: hmax, hmin
integer :: NewKZ, OldKZ, nlayers
!Begin-----------------------------------------------------------------
!WorkSize
ILB = Me%WorkSize%ILB
IUB = Me%WorkSize%IUB
JLB = Me%WorkSize%JLB
JUB = Me%WorkSize%JUB
KUB = Me%WorkSize%KUB
!Gets a pointer to the Bathymetry
call GetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructKFloor - Geometry - ERR01'
!Gets a pointer to 2D WaterPoints
call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructKFloor - Geometry - ERR02'
!Exterior Boundaries
call GetExteriorBoundaryFaces(Me%ObjHorizontalMap, &
BoundaryPointsFaceU = ExteriorFacesU, &
BoundaryPointsFaceV = ExteriorFacesV, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR03'
!Checks which is the domain close to the bottom
if (Me%FirstDomain%DomainType == FixSpacing .or. &
Me%FirstDomain%DomainType == FixSediment .or. &
Me%FirstDomain%DomainType == CartesianTop) then
do j = JLB, JUB
do i = ILB, IUB
if (WaterPoints2D(i, j) == WaterPoint) then
Me%KFloor%Domain(i, j) = Me%FirstDomain%ID
endif
enddo
enddo
else
do j = JLB, JUB
do i = ILB, IUB
if (WaterPoints2D(i, j) == WaterPoint) then
Me%KFloor%Domain(i, j) = Me%LastDomain%ID
CurrentDomain => Me%FirstDomain
do while (associated(CurrentDomain))
if (Bathymetry(i, j) > CurrentDomain%DomainDepth) then
Me%KFloor%Domain(i, j) = CurrentDomain%ID
exit
endif
CurrentDomain => CurrentDomain%Next
enddo
endif
enddo
enddo
endif
!Computes for KFloor%Z
NeedToStop = .false.
MaxDomainThickness = null_real
doj: do j = JLB, JUB
doi: do i = ILB, IUB
iw: if (WaterPoints2D(i, j) == WaterPoint) then
LayersBelow = 0
CurrentDomain => Me%FirstDomain
do while (CurrentDomain%ID /= Me%KFloor%Domain(i, j))
LayersBelow = LayersBelow + CurrentDomain%NumberOfLayers
CurrentDomain => CurrentDomain%Next
enddo
!Set Current domain so that it points to the domain close to the bottom.
CurrentDomain => Me%FirstDomain
do while (CurrentDomain%ID /= Me%KFloor%Domain(i, j))
CurrentDomain => CurrentDomain%Next
enddo
!In the case of a sigma or fixspacing domain
!KFloorZ is always equal to LayersBelow + 1
if ((CurrentDomain%DomainType /= Cartesian ) .and. &
! (CurrentDomain%DomainType /= Harmonic ) .and. &
(CurrentDomain%DomainType /= CartesianTop)) then
Me%KFloor%Z(i, j) = 1 + LayersBelow
else
!Domain Top Depth
if (CurrentDomain%ID == Me%LastDomain%ID) then
if (CurrentDomain%DomainType /= CartesianTop) then
TopDepth = 0.
else
if (present(SurfaceElevation)) then
TopDepth = dble(-1.0 * SurfaceElevation(i, j))
else
TopDepth = 0.
endif
endif
else
TopDepth = dble(CurrentDomain%DomainDepth)
endif
!Domain Bottom Depth
if (CurrentDomain%ID > 1) then
BottomDepth = dble(CurrentDomain%Prev%DomainDepth)
else
if (CurrentDomain%DomainType /= CartesianTop) then
call GetMaximumValue(Me%ObjTopography, Aux4)
BottomDepth = dble(Aux4)
else
!To minimize roundoff errors when BathymTopoFactor = 1
if (Me%ExternalVar%BathymTopoFactor/=1.0) then
BottomDepth = dble(Me%ExternalVar%BathymTopoFactor) * dble(Bathymetry(i, j))
else
BottomDepth = dble(Bathymetry(i, j))
endif
endif
endif
!Domain Thickness
DomainThickness = BottomDepth - TopDepth
!Verifies if the depth of the Bathymetry is between the upper and the lower
!part of the layer
FoundKFloor = .false.
LayerTopDepth = TopDepth
iLayer = CurrentDomain%UpperLayer
do while (iLayer >= CurrentDomain%LowerLayer .and. .not. FoundKFloor)
!MinimalThickness is now defined as a percentage of the initial thickness of each layer
MinimalThickness = dble(CurrentDomain%MinInitialLayerThickness)
!In the case of Cartesian, Harmonic domain
!the layer thickness is given in meters
if (CurrentDomain%DomainType == Cartesian .or. &
! CurrentDomain%DomainType == Harmonic .or. &
CurrentDomain%DomainType == CartesianTop) then
LayerThicknessMax = dble(CurrentDomain%LayerThickness(iLayer))
LayerThicknessMin = dble(CurrentDomain%LayerThickness(iLayer)) * &
MinimalThickness
else
LayerThicknessMax = dble(CurrentDomain%LayerThickness(iLayer)) * DomainThickness
LayerThicknessMin = dble(CurrentDomain%LayerThickness(iLayer)) * &
DomainThickness * MinimalThickness
endif
AllmostZero_ = AllmostZeroFraction * CurrentDomain%LayerThickness(iLayer)
LayerBottomDepthMin = LayerTopDepth + LayerThicknessMin - AllmostZero_
LayerBottomDepthMax = LayerTopDepth + LayerThicknessMax + AllmostZero_
AuxDepth = 0.0
!To minimize roundoff errors when BathymTopoFactor = 1
if (Me%ExternalVar%BathymTopoFactor/=1.0) then
AuxDepth = dble(Me%ExternalVar%BathymTopoFactor) * dble(Bathymetry(i, j))
else
AuxDepth = dble(Bathymetry(i, j))
endif
if (LayerBottomDepthMin <= AuxDepth .and. LayerBottomDepthMax >= AuxDepth) then
Me%KFloor%Z(i, j) = iLayer
FoundKFloor = .true.
!else if (Me%BathymNotCorrect .and. LayerBottomDepthMin > AuxDepth) then
else if (LayerBottomDepthMin > AuxDepth) then
Me%KFloor%Z(i, j) = iLayer + 1
FoundKFloor = .true.
else
LayerTopDepth = LayerTopDepth + LayerThicknessMax
endif
iLayer = iLayer - 1
enddo
!Points with Bathymetry > TopDepth
if (.not. FoundKFloor) then
if (CurrentDomain%DomainType /= CartesianTop) then
Me%KFloor%Z(i, j) = CurrentDomain%UpperLayer
else
MaxDomainThickness = max(MaxDomainThickness, DomainThickness)
TotalLayerThickness = LayerTopDepth - TopDepth
NeedToStop = .true.
endif
endif
endif
else iw
Me%KFloor%Z(i, j) = FillValueInt
endif iw
enddo doi
enddo doj
if (NeedToStop) then
write(*,*)'Domain Thickness Exceeds Layer Definition'
write(*,*)'Maximum DomainThickness : ',MaxDomainThickness
write(*,*)'Total Layer Thickness : ',TotalLayerThickness
stop 'ConstructKFloor - ModuleGeometry - ERR99'
endif
if (Me%LastDomain%DomainType == Sigma) then
if (Me%LastDomain%SigmaZleveHybrid) then
hmax = Me%LastDomain%SigmaZleveHybrid_Hmax
hmin = Me%LastDomain%SigmaZleveHybrid_Hmin
do j = JLB, JUB
do i = ILB, IUB
if (WaterPoints2D(i, j) == WaterPoint .and. Bathymetry(i,j) < hmax) then
OldKZ = Me%KFloor%Z(i, j)
nlayers = KUB - OldKZ
if (Bathymetry(i,j) > hmin) then
NewKZ = OldKZ + int(nlayers * (hmax - Bathymetry(i,j))/(hmax-hmin))
else
NewKZ = KUB
endif
if (NewKZ > OldKZ) then
Me%KFloor%Z(i, j) = NewKZ
endif
endif
enddo
enddo
endif
endif
!Computes KFloor%U
do j = JLB, JUB + 1
do i = ILB, IUB
!Water in (i,j) and in (i,j-1)? or current face is an exterior face
if ( ((Me%KFloor%Z(i, j ) > FillValueInt) .and. &
(Me%KFloor%Z(i, j - 1) > FillValueInt)) .or. &
(ExteriorFacesU(i,j) == 1) ) then
Me%KFloor%U(i, j) = max(Me%KFloor%Z(i, j), Me%KFloor%Z(i, j - 1))
else
Me%KFloor%U(i, j) = FillValueInt
endif
enddo
enddo
!Computes KFloor%V
do j = JLB, JUB
do i = ILB, IUB + 1
!Water in (i,j) and in (i-1,j)? or current face is an exterior face
if ( ((Me%KFloor%Z(i, j) > FillValueInt) .and. &
(Me%KFloor%Z(i - 1, j) > FillValueInt)) .or. &
(ExteriorFacesV(i,j) == 1)) then
Me%KFloor%V(i, j) = max(Me%KFloor%Z(i, j), Me%KFloor%Z(i - 1, j))
else
Me%KFloor%V(i, j) = FillValueInt
endif
enddo
enddo
call UnGetHorizontalMap(Me%ObjHorizontalMap, ExteriorFacesU, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR04'
call UnGetHorizontalMap(Me%ObjHorizontalMap, ExteriorFacesV, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR05'
!Disposes pointer to the Bathymetry
call UngetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR06'
!Disposes pointer to WaterPoints2D
call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
if (STAT_CALL .NE. SUCCESS_) stop 'ConstructKFloor - ModuleGeometry - ERR07'
if (Me%RemoveLandBotLayersON) then
call RemoveLandBottomLayers
endif
!----------------------------------------------------------------------
end subroutine ConstructKFloor
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine Reallocate3DGeometryPropR8(GeoProp3D, Matrix3D, KLB, KUB, KLB_Max, KUB_Max)
!Arguments-------------------------------------------------------------
real(8), allocatable, dimension(:,:,:) :: GeoProp3D
real(8), pointer, dimension(:,:,:) :: Matrix3D
integer :: KLB, KUB, KLB_Max, KUB_Max
!Local-----------------------------------------------------------------
integer :: STATUS
integer :: ILB, IUB, JLB, JUB
!Begin-----------------------------------------------------------------
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
Matrix3D(:,:,:) = GeoProp3D(:,:,:)
deallocate (GeoProp3D, stat = STATUS)
if (STATUS /= SUCCESS_) stop 'Reallocate3DGeometryPropR8 - Geometry - ERR10'
allocate (GeoProp3D(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'Reallocate3DGeometryPropR8 - Geometry - ERR20'
GeoProp3D(:,:,:) = Matrix3D(:,:,KLB_Max:KUB_Max)
end subroutine Reallocate3DGeometryPropR8
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine Reallocate3DGeometryProp(GeoProp3D, Matrix3D, KLB, KUB, KLB_Max, KUB_Max)
!Arguments-------------------------------------------------------------
real, allocatable, dimension(:,:,:) :: GeoProp3D
real(8), pointer, dimension(:,:,:) :: Matrix3D
integer :: KLB, KUB, KLB_Max, KUB_Max
!Local-----------------------------------------------------------------
integer :: STATUS
integer :: ILB, IUB, JLB, JUB
!Begin-----------------------------------------------------------------
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
Matrix3D(:,:,:) = GeoProp3D(:,:,:)
deallocate (GeoProp3D, stat = STATUS)
if (STATUS /= SUCCESS_) stop 'Reallocate3DGeometryProp - Geometry - ERR10'
allocate (GeoProp3D(ILB:IUB, JLB:JUB, KLB:KUB), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'Reallocate3DGeometryProp - Geometry - ERR20'
GeoProp3D(:,:,:) = Matrix3D(:,:,KLB_Max:KUB_Max)
end subroutine Reallocate3DGeometryProp
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine Reallocate1DGeometryProp(GeoProp1D, KLD, KUD, KLD_New, KUD_New)
!Arguments-------------------------------------------------------------
real, pointer, dimension(:) :: GeoProp1D
integer :: KLD, KUD, KLD_New, KUD_New
!Local-----------------------------------------------------------------
real, pointer, dimension(:) :: Matrix1D
integer :: STATUS
!Begin-----------------------------------------------------------------
allocate(Matrix1D(KLD:KUD))
Matrix1D(:) = GeoProp1D(:)
deallocate (GeoProp1D, stat = STATUS)
if (STATUS /= SUCCESS_) stop 'Reallocate1DGeometryProp - Geometry - ERR10'
allocate (GeoProp1D(KLD_New:KUD_New), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'Reallocate1DGeometryProp - Geometry - ERR20'
GeoProp1D(KLD_New:KUD_New) = Matrix1D(KLD:KUD)
deallocate(Matrix1D)
end subroutine Reallocate1DGeometryProp
!--------------------------------------------------------------------------
subroutine RemoveLandBottomLayers
!Parameter-------------------------------------------------------------
!Local-----------------------------------------------------------------
real(8), pointer, dimension(:,:,:) :: Matrix3D
integer :: STATUS
integer :: ILB, IUB, JLB, JUB, KLB, KUB
integer :: i, j, Kmin, KLB_Max, KUB_Max
type (T_Domain), pointer :: CurrentDomain
integer :: KLD_New, KUD_New, KLD, KUD, KminAux
!Begin-----------------------------------------------------------------
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
!Find Kmin
Kmin = - FillValueInt
do j = JLB, JUB
do i = ILB, IUB
if (Me%KFloor%Z(i, j) > 0) then
if (Me%KFloor%Z(i, j) < Kmin) then
Kmin = Me%KFloor%Z(i, j)
endif
endif
enddo
enddo
KminAux = Kmin
#ifdef _USE_MPI
call GetKfloorZminMPI(HorizontalGridID = Me%ObjHorizontalGrid, &
KminZin = KminAux, &
KminZout = Kmin, &
STAT = STATUS)
if (STATUS /= SUCCESS_) stop 'RemoveLandBottomLayers - Geometry - ERR10'
#endif
!Correct the 2D Variables
do j = JLB, JUB
do i = ILB, IUB
if (Me%KFloor%Z(i, j) > 0) then
Me%KFloor%Z(i, j) = Me%KFloor%Z(i, j) - Kmin + 1
endif
enddo
enddo
do j = JLB, JUB
do i = ILB, IUB
if (Me%KFloor%U(i, j) > 0) then
Me%KFloor%U(i, j) = Me%KFloor%U(i, j) - Kmin + 1
endif
enddo
enddo
do j = JLB, JUB
do i = ILB, IUB
if (Me%KFloor%V(i, j) > 0) then
Me%KFloor%V(i, j) = Me%KFloor%V(i, j) - Kmin + 1
endif
enddo
enddo
!reallocate and shift vertically the 3D Variables
KLB = Me%Size%KLB
KUB_Max = Me%Size%KUB
allocate(Matrix3D(ILB:IUB, JLB:JUB, KLB:KUB_Max), stat = STATUS)
if (STATUS /= SUCCESS_) stop 'RemoveLandBottomLayers - Geometry - ERR20'
Me%Size%KUB = Me%Size%KUB - Kmin + 1
Me%WorkSize%KUB = Me%WorkSize%KUB - Kmin + 1
KUB = Me%Size%KUB
KLB_Max = KLB + Kmin -1
!VolumeZ
call Reallocate3DGeometryPropR8(GeoProp3D = Me%Volumes%VolumeZ, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!VolumeU
call Reallocate3DGeometryPropR8(GeoProp3D = Me%Volumes%VolumeU, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!VolumeV
call Reallocate3DGeometryPropR8(GeoProp3D = Me%Volumes%VolumeV, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!VolumeZOld
call Reallocate3DGeometryPropR8(GeoProp3D = Me%Volumes%VolumeZOld, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!AreaU
call Reallocate3DGeometryProp(GeoProp3D = Me%Areas%AreaU, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!AreaV
call Reallocate3DGeometryProp(GeoProp3D = Me%Areas%AreaV, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!SZZ
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%SZZ, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DZZ
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DZZ, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DWZ
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DWZ, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DUZ
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DUZ, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DVZ
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DVZ, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!InitialSZZ
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%InitialSZZ, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!ZCellCenter
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%ZCellCenter, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DZI
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DZI, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DZE
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DZE, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DWZ_Xgrad
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DWZ_Xgrad, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!DWZ_Ygrad
call Reallocate3DGeometryProp(GeoProp3D = Me%Distances%DWZ_Ygrad, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
!NearbyAvgVel_Z
call Reallocate3DGeometryProp(GeoProp3D = Me%NearbyAvgVel_Z, &
Matrix3D = Matrix3D, &
KLB = KLB , &
KUB = KUB , &
KLB_Max = KLB_Max , &
KUB_Max = KUB_Max )
CurrentDomain => Me%FirstDomain
do while (associated(CurrentDomain))
KLD = CurrentDomain%LowerLayer
KUD = CurrentDomain%UpperLayer
KLD_New = KLD - Kmin + 1
KUD_New = KUD - Kmin + 1
call Reallocate1DGeometryProp(GeoProp1D = CurrentDomain%LayerThickness, &
KLD = KLD , &
KUD = KUD , &
KLD_New = KLD_New, &
KUD_New = KUD_New)
call Reallocate1DGeometryProp(GeoProp1D = CurrentDomain%LayerMinThickness, &
KLD = KLD , &
KUD = KUD , &
KLD_New = KLD_New, &
KUD_New = KUD_New)
call Reallocate1DGeometryProp(GeoProp1D = CurrentDomain%LayerMaxThickness, &
KLD = KLD , &
KUD = KUD , &
KLD_New = KLD_New, &
KUD_New = KUD_New)
!Atualizes UpperLayer and LowerLayer of each domain
CurrentDomain%LowerLayer = KLD_New
CurrentDomain%UpperLayer = KUD_New
CurrentDomain => CurrentDomain%Next
enddo
deallocate(Matrix3D, stat = STATUS)
if (STATUS /= SUCCESS_) stop 'RemoveLandBottomLayers - Geometry - ERR30'
end subroutine RemoveLandBottomLayers
!--------------------------------------------------------------------------
subroutine Add_Domain(NewDomain)
!Arguments-------------------------------------------------------------
type(T_Domain), pointer :: NewDomain
!----------------------------------------------------------------------
! Add to the WaterDomain List a new property
if (.not.associated(Me%FirstDomain)) then
Me%FirstDomain => NewDomain
Me%LastDomain => NewDomain
else
NewDomain%Prev => Me%LastDomain
Me%LastDomain%Next => NewDomain
Me%LastDomain => NewDomain
end if
!----------------------------------------------------------------------
end subroutine Add_Domain
!--------------------------------------------------------------------------
#ifdef _USE_SEQASSIMILATION
subroutine PointToGeometryState(GeometryID, STAT)
!Arguments---------------------------------------------------------------
integer :: GeometryID
integer, optional, intent(OUT) :: STAT
!Local-------------------------------------------------------------------
integer :: ready_
integer :: STAT_
!------------------------------------------------------------------------
!This is a compilation of points (one for each variable) for internal memory spaces
STAT_ = UNKNOWN_
call Ready(GeometryID, ready_)
cd1 : if (ready_ .EQ. IDLE_ERR_) then
Me%AuxPointer%SZZ => Me%Distances%SZZ
Me%AuxPointer%DWZ => Me%Distances%DWZ
Me%AuxPointer%DUZ => Me%Distances%DUZ
Me%AuxPointer%DVZ => Me%Distances%DVZ
Me%AuxPointer%DZZ => Me%Distances%DZZ
Me%AuxPointer%ZCellCenter => Me%Distances%ZCellCenter
Me%AuxPointer%AreaU => Me%Areas%AreaU
Me%AuxPointer%AreaV => Me%Areas%AreaV
Me%AuxPointer%WaterColumnU => Me%WaterColumn%U
Me%AuxPointer%WaterColumnV => Me%WaterColumn%V
Me%AuxPointer%WaterColumnZ => Me%WaterColumn%Z
Me%AuxPointer%VolumeZ => Me%Volumes%VolumeZ
Me%AuxPointer%VolumeU => Me%Volumes%VolumeU
Me%AuxPointer%VolumeV => Me%Volumes%VolumeV
Me%AuxPointer%VolumeZOld => Me%Volumes%VolumeZOld
STAT_ = SUCCESS_
else
STAT_ = ready_
end if cd1
if (present(STAT))STAT = STAT_
end subroutine PointToGeometryState
#endif _USE_SEQASSIMILATION
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MODIFIER MO
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine ComputeInitialGeometry(GeometryID, WaterPoints3D, SurfaceElevation, &
ActualTime, ContinuesCompute, NonHydrostatic, &
SZZ, STAT)
!Arguments-------------------------------------------------------------
integer :: GeometryID
integer, dimension(:, :, :), pointer :: WaterPoints3D
real, dimension(:, :), pointer, optional :: SurfaceElevation
type(T_Time), optional :: ActualTime
logical, intent(in), optional :: ContinuesCompute
logical, intent(in), optional :: NonHydrostatic
real, dimension(:, :, :), pointer, optional :: SZZ
integer, intent(out), optional :: STAT
!Local-----------------------------------------------------------------
integer :: ready_
integer :: STAT_
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
call Ready(GeometryID, ready_)
cd1 : if (ready_ .EQ. IDLE_ERR_) then
!Actualize the time
if (present(ActualTime)) then
Me%ActualTime = ActualTime
endif
if (present(NonHydrostatic)) then
Me%ExternalVar%NonHydrostatic = NonHydrostatic
else
Me%ExternalVar%NonHydrostatic = .false.
endif
if (present(ContinuesCompute)) then
Me%ExternalVar%ContinuesCompute = ContinuesCompute
else
Me%ExternalVar%ContinuesCompute = .false.
endif
cd2 : if (Me%ExternalVar%ContinuesCompute) then
!Computes the Distances
call ComputeDistances (WaterPoints3D)
!Computes the Areas
call ComputeAreas
!Computes the Volumes
call ComputeVolumes
!It is necessary for the Soil model
call ComputeZCellCenter
else cd2
if (present(SZZ)) then
call SetMatrixValue( GetPointer(Me%Distances%SZZ), Me%Size, SZZ )
else
!Constructs SZZ with the initial surface elevation
call ComputeSZZ (SurfaceElevation, INITIALGEOMETRY, WaterPoints3D = WaterPoints3D)
if (Me%LastDomain%DomainType == Sigma) then
if (Me%LastDomain%SigmaZleveHybrid) then
! SigmaZHybridSurface
call SigmaZHybridSurface(SurfaceElevation, WaterPoints3D)
endif
endif
endif
!Computes the Distances
call ComputeDistances (WaterPoints3D)
!Computes the Volumes
call ComputeVolumes
!Computes the Areas
call ComputeAreas
!Stores VolumeZOLD
!In this case VolumeZOld will be equal to VolumeZ
call StoreVolumeZOld
!It is necessary for the Soil model
call ComputeZCellCenter
end if cd2
!Computes the WaterColumn
if (present(SurfaceElevation)) then
call ComputeWaterColumn (SurfaceElevation)
call ComputeVolume2D
endif
STAT_ = SUCCESS_
else cd1
STAT_ = ready_
end if cd1
if (present(STAT)) STAT = STAT_
end subroutine ComputeInitialGeometry
!--------------------------------------------------------------------------
subroutine SigmaZHybridSurface (SurfaceElevation, WaterPoints3D)
!Arguments-------------------------------------------------------------
real, dimension(:,: ), pointer :: SurfaceElevation
integer, dimension(:,:,:), pointer :: WaterPoints3D
!Local-----------------------------------------------------------------
integer :: STAT_CALL
integer :: ILB, IUB, JLB, JUB, KLB, KUB, kbottom
integer :: i, j, k
real(8) :: hmax
real(8) :: TotalThickness, TotalWaterColumn, AuxRacio
real, dimension(:,:), pointer :: Bathymetry
!Begin-----------------------------------------------------------------
ILB = Me%WorkSize%ILB
IUB = Me%WorkSize%IUB
JLB = Me%WorkSize%JLB
JUB = Me%WorkSize%JUB
KLB = Me%WorkSize%KLB
KUB = Me%WorkSize%KUB
!Gets Bathymetry
call GetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'SigmaZHybridSurface - Geometry - ERR10'
hmax = Me%LastDomain%SigmaZleveHybrid_Hmax
do j = JLB, JUB
do i = ILB, IUB
if (WaterPoints3D(i, j, KUB) == WaterPoint .and. Bathymetry(i,j) < hmax) then
kbottom = Me%KFloor%Z(i, j)
TotalThickness = sum(Me%LastDomain%LayerThickness(kbottom:KUB))
TotalWaterColumn = Bathymetry(i, j) + SurfaceElevation(i, j)
AuxRacio = TotalWaterColumn / TotalThickness
do k = KUB,KLB-1,-1
!surface
if (k == KUB) then
Me%Distances%SZZ(i, j, k) = -1.* SurfaceElevation(i, j)
!water column
else if (k >=kbottom .and. k < KUB) then
Me%Distances%SZZ(i, j, k) = Me%Distances%SZZ(i, j, k+1) + Me%LastDomain%LayerThickness(k+1) * AuxRacio
!bottom
else if (k == kbottom - 1) then
Me%Distances%SZZ(i, j, k) = Bathymetry(i, j)
!land
else
Me%Distances%SZZ(i, j, k) = FillValueReal
endif
enddo
endif
enddo
enddo
!UnGets Bathymetry
call UnGetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'SigmaZHybridSurface - Geometry - ERR20'
end subroutine SigmaZHybridSurface
!-------------------------------------------------------------------------------
subroutine UpdateKfloor(GeometryID, SurfaceElevation, BathymNotCorrect, STAT)
!Arguments-------------------------------------------------------------
integer :: GeometryID
real, dimension(:, :), pointer, optional :: SurfaceElevation
logical, optional :: BathymNotCorrect
integer, intent(out), optional :: STAT
!Local-----------------------------------------------------------------
integer :: ready_
integer :: STAT_
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
call Ready(GeometryID, ready_)
cd1 : if (ready_ .EQ. IDLE_ERR_) then
if (present(BathymNotCorrect)) then
Me%BathymNotCorrect = BathymNotCorrect
else
Me%BathymNotCorrect = .false.
endif
!Updates the Matrixes which contains the Ks - KFloor, etc...
if (present(SurfaceElevation)) then
call ConstructKFloor (SurfaceElevation)
else
call ConstructKFloor
endif
STAT_ = SUCCESS_
else cd1
STAT_ = ready_
end if cd1
if (present(STAT)) STAT = STAT_
end subroutine UpdateKfloor
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! This subroutine computes the VerticalGeometry
! It should be called during the working cycle
subroutine ComputeVerticalGeometry(GeometryID, WaterPoints3D, SurfaceElevation, &
ActualTime, VerticalVelocity, DT_Waterlevel, &
SZZ, DecayTime, KTop, OpenPoints3D, STAT)
!Arguments-------------------------------------------------------------
integer :: GeometryID
integer, dimension(:, :, :), pointer :: WaterPoints3D
integer, dimension(:, :, :), pointer, optional :: OpenPoints3D
real, dimension(:, :), pointer, optional :: SurfaceElevation
type (T_Time), optional :: ActualTime
real, dimension(:, :, :), pointer, optional :: VerticalVelocity !Gives the vertical variation
real, intent(in), optional :: DT_Waterlevel !for the lagragean coordinate
real, dimension(:, :, :), pointer, optional :: SZZ, DecayTime
integer, dimension(:, :), pointer, optional :: KTop
integer, intent(out), optional :: STAT
!Local-----------------------------------------------------------------
integer :: ready_
integer :: STAT_
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
call Ready(GeometryID, ready_)
cd1 : if (ready_ .EQ. IDLE_ERR_) then
!Actualize the time
if (present(ActualTime)) then
Me%ActualTime = ActualTime
endif
if (present(SZZ)) then
call SetMatrixValue( GetPointer(Me%Distances%SZZ), Me%Size, SZZ )
else
!Computes SZZ
call ComputeSZZ(SurfaceElevation, TRANSIENTGEOMETRY, VerticalVelocity, DT_Waterlevel, WaterPoints3D, &
OpenPoints3D)
if (Me%LastDomain%DomainType == Sigma) then
if (Me%LastDomain%SigmaZleveHybrid) then
! SigmaZHybridSurface
call SigmaZHybridSurface(SurfaceElevation, WaterPoints3D)
endif
endif
endif
if(present(KTop))then
Me%KTop%Z(:,:) = KTop(:,:)
endif
if(present(DecayTime)) then
Me%ExternalVar%DecayTime => DecayTime
endif
!Stores VolumeZOld
call StoreVolumeZOld
!Computes Distances
call ComputeDistances(WaterPoints3D)
!Computes Volumes
call ComputeVolumes
!Computes Areas
call ComputeAreas
!It is necessary for the Soil model
call ComputeZCellCenter
if (present(SurfaceElevation)) then
!Computes the WaterColumn
if (Me%FirstDomain%DomainType /= FixSediment) then
call ComputeWaterColumn(SurfaceElevation)
endif
endif
nullify(Me%Externalvar%DecayTime)
STAT_ = SUCCESS_
else cd1
STAT_ = ready_
end if cd1
if (present(STAT)) STAT = STAT_
end subroutine ComputeVerticalGeometry
!--------------------------------------------------------------------------
subroutine ComputeKTop (SedimentDomain)
!Arguments-------------------------------------------------------------
type(T_Domain), pointer :: SedimentDomain
!External--------------------------------------------------------------
integer :: STAT_CALL
!Local-----------------------------------------------------------------
integer, dimension(:,:), pointer :: WaterPoints2D
integer :: ILB, IUB, JLB, JUB
integer :: i, j
!Begin-----------------------------------------------------------------
ILB = Me%WorkSize%ILB
IUB = Me%WorkSize%IUB
JLB = Me%WorkSize%JLB
JUB = Me%WorkSize%JUB
!Gets a pointer to 2D WaterPoints
call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ComputeKTop - ModuleGeometry - ERR01'
do j = JLB, JUB
do i = ILB, IUB
if (WaterPoints2D(i, j) == WaterPoint) then
Me%KTop%Z(i, j) = Me%WorkSize%KUB - SedimentDomain%EmptyTopLayers
else
Me%KTop%Z(i, j) = Me%WorkSize%KUB
endif
enddo
enddo
!Disposes pointer to WaterPoints2D
call UnGetHorizontalMap(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ComputeKTop - ModuleGeometry - ERR02'
end subroutine ComputeKTop
!--------------------------------------------------------------------------
!Computes the WaterColumn (Bathymetry + SurfaceElevation)
subroutine ComputeWaterColumn(SurfaceElevation)
!Arguments-------------------------------------------------------------
real, dimension(:, :), pointer :: SurfaceElevation
!Local-----------------------------------------------------------------
real, dimension(:,:,:), pointer :: DUZ, DVZ
real, dimension(:, : ), pointer :: Bathymetry, WaterColumnZ, &
WaterColumnU, WaterColumnV
integer, dimension(:, :), pointer :: WaterPoints2D, KFloorU, KFloorV
integer :: ILB, IUB, JLB, JUB, KUB
integer :: i, j, k, STAT_CALL, kbottom
integer :: CHUNK
!Worksize
ILB = Me%WorkSize%ILB
IUB = Me%WorkSize%IUB
JLB = Me%WorkSize%JLB
JUB = Me%WorkSize%JUB
KUB = Me%WorkSize%KUB
WaterColumnZ => Me%WaterColumn%Z
WaterColumnU => Me%WaterColumn%U
WaterColumnV => Me%WaterColumn%V
KFloorU => Me%KFloor%U
KFloorV => Me%KFloor%V
DUZ => Me%Distances%DUZ
DVZ => Me%Distances%DVZ
!Gets Bathymetry
call GetGridData(Me%ObjTopography, Bathymetry, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ComputeWaterColumn - Geometry - ERR01'
!Gets WaterPoints2D
call GetWaterPoints2D(Me%ObjHorizontalMap, WaterPoints2D, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
stop 'ComputeWaterColumn - Geometry - ERR02'
if (MonitorPerformance) call StartWatch ("ModuleGeometry", "ComputeWaterColumn")
CHUNK = Chunk_J(JLB,JUB)
!Computes WaterColumn
!$OMP PARALLEL PRIVATE(i,j,k,kbottom)
!$OMP DO SCHEDULE(DYNAMIC, CHUNK)
do j = JLB, JUB
do i = ILB, IUB
if (WaterPoints2D(i, j) == WaterPoint) &
WaterColumnZ(i, j)