Skip to content
Permalink
Branch: master
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
18300 lines (12724 sloc) 685 KB
!------------------------------------------------------------------------------
! IST/MARETEC, Water Modelling Group, Mohid modelling system
!------------------------------------------------------------------------------
!
! TITLE : Mohid Model
! PROJECT : Mohid Base 2
! MODULE : HorizontalGrid
! URL : http://www.mohid.com
! AFFILIATION : IST/MARETEC, Marine Modelling Group
! DATE : May 2003
! REVISION : Frank Braunschweig - v4.0
! DESCRIPTION : Module to read and calculate a horizontal Grid
!
!------------------------------------------------------------------------------
!
!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 ModuleHorizontalGrid
use ModuleGlobalData
use ModuleFunctions, only : Rodaxy, FromCartesianToGrid, FromGridToCartesian, &
UTMToLatLon, LatLonToUTM, ComputeGridZone, &
LatLonToLambertSP2, RelativePosition4VertPolygon, &
CHUNK_J, WGS84toGoogleMaps, AngleFromFieldToGrid, &
AngleFromGridToField, THOMAS_2D, THOMAS_3D
#ifdef _USE_PROJ4
use ModuleFunctions, only : GeographicToCartesian, CartesianToGeographic
#endif _USE_PROJ4
#ifdef _USE_MPI
use ModuleFunctions, only : MPIKind
use mpi
#endif _USE_MPI
use ModuleStopWatch, only : StartWatch, StopWatch
use ModuleEnterData
use ModuleDrawing
use ModuleHDF5
implicit none
private
!Subroutine----------------------------------------------------------------
!Constructor
public :: ConstructHorizontalGrid
private :: AllocateInstance
private :: ConstructGlobalVariables
private :: AllocateVariables
private :: Mercator
private :: ConvertCoordenates !Fromer CONVCO_BAT
private :: GridPointsLocation1D
private :: GridPointsLocation2D
private :: ComputeDistances
private :: ComputeRotation
private :: ComputeGridCellArea
private :: ComputeCoriolis
public :: ConstructFatherGridLocation
private :: ConstructNewFatherGrid1D
private :: ConstructNewFatherGrid2D
private :: ConstructIWDVel
private :: DetermineMaxRatio
private :: Add_FatherGrid
private :: CheckGridBorder
private :: DefineBorderPolygons
private :: DefinesBorderPoly
public :: ConstructIWDTwoWay
!Modifier
public :: WriteHorizontalGrid
public :: WriteHorizontalGrid_UV
public :: LocateCell
public :: LocateCell1D
public :: LocateCellPolygons
public :: RecenterHorizontalGrid
public :: Add_MPI_ID_2_Filename
public :: No_BoxMap_HaloArea
!Selector
public :: GetHorizontalGridSize
public :: GetHorizontalGrid
public :: GetGridOrigin
public :: GetGridAngle
public :: GetGridZone
public :: GetLatitudeLongitude
public :: GetGridCoordType
public :: GetCoordTypeList
public :: GetGridMeanLatLong
public :: GetCoriolisFrequency
public :: GetGridLatitudeLongitude
public :: GetComputeZUV
public :: GetGridFileName
public :: GetCheckDistortion
public :: GetGridRotation
public :: GetCellRotation
public :: GetDefineCellsMap
public :: GetNotDefinedCells
public :: GetZCoordinates
public :: GetCornersCoordinates
public :: GetFatherGridID
public :: GetGridCellArea
private :: Search_FatherGrid
public :: GetXYCellZ
public :: GetXYCellZ_ThreadSafe
public :: GetCellZ_XY
public :: GetCellZInterceptByLine
public :: GetCellZInterceptByPolygon
public :: GetCellZInterceptByXYZPoints
public :: GetGridBorderPolygon
public :: GetGridOutBorderPolygon
public :: GetGridBorderCartPolygon
public :: GetGridOutBorderCartPolygon
public :: GetGridBorderLimits
public :: GetGridOutBorderCartLimits
public :: GetXYInsideDomain
private :: InsideDomainPolygon
public :: GetGridBorderType
public :: GetDDecompOpenBorders
public :: GetDDecompParameters
public :: GetDDecompSlaves
public :: GetDDecompSlavesSize
public :: GetDDecompWorkSize2D
public :: GetDDecompMapping2D
public :: GetDDecompMPI_ID
public :: GetDDecompON
public :: GetSonWindow
public :: GetTwoWayAux
public :: UnGetTwoWayAux
public :: UnGetHorizontalGrid
public :: InterpolXYPoint
public :: ReturnsIntersectionCorners
private :: ComputesIntersectionCorners
public :: WindowIntersectDomain
private :: WindowCellsIntersection
#ifdef _USE_PROJ4
public :: FromGeo2SpherMercator1D
public :: FromGeo2SpherMercatorScalar
public :: FromGeographic2SphericMercator
#endif
!Destructor
public :: KillHorizontalGrid
private :: DeallocateInstance
private :: KillFatherGridList
!Management
private :: Ready
!Auxilary
private :: USCONVCO
private :: UTGP83
private :: GPUT83
private :: IUTPC
private :: IUTGP
private :: DRGPUT
private :: DRUTGP
private :: TODMS
private :: DATUMM
private :: TMGRID
private :: TCONST
private :: TCONPC
private :: TMGEOD
private :: wgs84_to_rd
private :: wgs842bessel_
private :: bessel2rd_
private :: rd_to_wgs84
private :: rd2bessel_
private :: bessel2wgs84_
!Interfaces----------------------------------------------------------------
interface ConstructHorizontalGrid
module procedure ConstructHorizontalGridV1
module procedure ConstructHorizontalGridV2
module procedure ConstructHorizontalGridV3
end interface ConstructHorizontalGrid
interface LocateCell1D
module procedure LocateCell1DR4
module procedure LocateCell1DR8
end interface LocateCell1D
private :: UnGetHorizontalGrid1d
private :: UnGetHorizontalGrid2d
private :: UngetHorizontalGrid1DInt
private :: UngetHorizontalGrid2DInt
private :: UnGetHorizontalGridPolygon
interface UnGetHorizontalGrid
module procedure UnGetHorizontalGrid1D
module procedure UnGetHorizontalGrid2D
module procedure UnGetHorizontalGrid1DInt
module procedure UnGetHorizontalGrid2DInt
module procedure UnGetHorizontalGridPolygon
end interface UnGetHorizontalGrid
public :: InterpolRegularGrid
interface InterpolRegularGrid
module procedure InterpolRegularGrid2D
module procedure InterpolRegularGrid3D
module procedure InterpolRegularGrid3D8
end interface InterpolRegularGrid
private:: RotateVectorFieldToGrid2D
private:: RotateVectorFieldToGrid3D
public :: RotateVectorFieldToGrid
interface RotateVectorFieldToGrid
module procedure RotateVectorFieldToGrid2D
module procedure RotateVectorFieldToGrid3D
end interface RotateVectorFieldToGrid
private:: RotateVectorGridToField2DR4
private:: RotateVectorGridToField2DR8
private:: RotateVectorGridToField3D
public :: RotateVectorGridToField
interface RotateVectorGridToField
module procedure RotateVectorGridToField2DR4
module procedure RotateVectorGridToField2DR8
module procedure RotateVectorGridToField3D
end interface RotateVectorGridToField
private:: RotateAngleFieldToGrid2D
private:: RotateAngleFieldToGrid3D
public :: RotateAngleFieldToGrid
interface RotateAngleFieldToGrid
module procedure RotateAngleFieldToGrid2D
module procedure RotateAngleFieldToGrid3D
end interface RotateAngleFieldToGrid
private:: RotateAngleGridToField2D
private:: RotateAngleGridToField3D
public :: RotateAngleGridToField
interface RotateAngleGridToField
module procedure RotateAngleGridToField2D
module procedure RotateAngleGridToField3D
end interface RotateAngleGridToField
private:: ComputeAngleFromGridComponents2D
private:: ComputeAngleFromGridComponents3D
public :: ComputeAngleFromGridComponents
interface ComputeAngleFromGridComponents
module procedure ComputeAngleFromGridComponents2D
module procedure ComputeAngleFromGridComponents3D
end interface ComputeAngleFromGridComponents
#ifdef _USE_MPI
public :: THOMAS_DDecompHorizGrid
interface THOMAS_DDecompHorizGrid
module procedure THOMAS_2D_DDecompHorizGrid
module procedure THOMAS_3D_DDecompHorizGrid
end interface
public :: JoinGridData
interface JoinGridData
module procedure JoinGridData_1D
module procedure JoinGridData_2D_R4
module procedure JoinGridData_2D_R8
module procedure JoinGridData_2Dint
module procedure JoinGridData_3D_R4
module procedure JoinGridData_3D_R8
module procedure JoinGridData_3Dint
end interface
public :: JoinGridData_In
interface JoinGridData_In
module procedure JoinGridData_1D_In
module procedure JoinGridData_2D_R4_In
module procedure JoinGridData_2D_R8_In
module procedure JoinGridData_2Dint_In
module procedure JoinGridData_3D_R4_In
module procedure JoinGridData_3D_R8_In
module procedure JoinGridData_3Dint_In
end interface
public :: BroadcastGridData
interface BroadcastGridData
module procedure BroadcastGridData2D_R4
module procedure BroadcastGridData2D_R8
module procedure BroadcastGridData3D_R4
module procedure BroadcastGridData3D_R8
end interface
public :: BroadcastGridData_In
interface BroadcastGridData_In
module procedure BroadcastGridData2D_R4_In
module procedure BroadcastGridData2D_R8_In
module procedure BroadcastGridData3D_R4_In
module procedure BroadcastGridData3D_R8_In
end interface
public :: ReceiveSendProperitiesMPI
interface ReceiveSendProperitiesMPI
module procedure ReceiveSendProperities3DMPIr4
module procedure ReceiveSendProperities2DMPIr4
module procedure ReceiveSendProperities3DMPIr8
module procedure ReceiveSendProperities2DMPIr8
endinterface
public :: ReceiveSendLogicalMPI
interface ReceiveSendLogicalMPI
module procedure AtLeastOneDomainIsTrue
endinterface
public :: GetKfloorZminMPI
#endif _USE_MPI
!Parameter-----------------------------------------------------------------
!Block information
character(LEN = StringLength), parameter :: BeginXX = '<BeginXX>'
character(LEN = StringLength), parameter :: EndXX = '<EndXX>'
character(LEN = StringLength), parameter :: BeginYY = '<BeginYY>'
character(LEN = StringLength), parameter :: EndYY = '<EndYY>'
character(LEN = StringLength), parameter :: BeginCornersXY = '<CornersXY>'
character(LEN = StringLength), parameter :: EndCornersXY = '<'//backslash//'CornersXY>'
character(LEN = StringLength), parameter :: BeginCartCornersXY = '<CartCornersXY>'
character(LEN = StringLength), parameter :: EndCartCornersXY = '<'//backslash//'CartCornersXY>'
!Calculation Points
integer, parameter :: ComputeZ_ = 1
integer, parameter :: ComputeU_ = 2
integer, parameter :: ComputeV_ = 3
integer, parameter :: ComputeCross_ = 4
integer, parameter :: ComputeZU_ = 5
integer, parameter :: ComputeZV_ = 6
integer, dimension(1:4), parameter :: diFace = (/0,1 ,0,-1/), djFace = (/-1, 0, 1, 0/)
integer, parameter :: LargeScaleModel_ = 1
integer, parameter :: Assimila_ = 2
integer, parameter :: SurfaceMM5_ = 3
!Input / Output
integer, parameter :: FileOpen = 1, FileClose = 0
!Type----------------------------------------------------------------------
type T_BorderLimits
real, dimension(4) :: Values = FillValueReal
logical :: ON = .false.
end type T_BorderLimits
type T_Compute
real, dimension(:), pointer :: XX_Z => null()
real, dimension(:), pointer :: YY_Z => null()
real, dimension(:), pointer :: XX_U => null()
real, dimension(:), pointer :: YY_U => null()
real, dimension(:), pointer :: XX_V => null()
real, dimension(:), pointer :: YY_V => null()
real, dimension(:), pointer :: XX_Cross => null()
real, dimension(:), pointer :: YY_Cross => null()
real, dimension(:,:), pointer :: XX2D_Z => null()
real, dimension(:,:), pointer :: YY2D_Z => null()
real, dimension(:,:), pointer :: XX2D_U => null()
real, dimension(:,:), pointer :: YY2D_U => null()
real, dimension(:,:), pointer :: XX2D_V => null()
real, dimension(:,:), pointer :: YY2D_V => null()
end type T_Compute
type T_Window
logical :: ON = .false. !initialization: jauch
integer :: ILB = null_int !initialization: jauch
integer :: IUB = null_int !initialization: jauch
integer :: JLB = null_int !initialization: jauch
integer :: JUB = null_int !initialization: jauch
end type T_Window
type T_FatherGrid
integer :: GridID = null_int !initialization: jauch - or should be set to 0 (zero)?
real, dimension(:,:), pointer :: XX_Z => null()
real, dimension(:,:), pointer :: YY_Z => null()
real, dimension(:,:), pointer :: XX_U => null()
real, dimension(:,:), pointer :: YY_U => null()
real, dimension(:,:), pointer :: XX_V => null()
real, dimension(:,:), pointer :: YY_V => null()
real, dimension(:,:), pointer :: XX_Cross => null()
real, dimension(:,:), pointer :: YY_Cross => null()
logical :: OkZ = .false. !initialization: jauch
logical :: OkU = .false. !initialization: jauch
logical :: OkV = .false. !initialization: jauch
logical :: OkCross = .false. !initialization: jauch
logical :: CornersXYInput = .false.
type (T_Window) :: Window
real :: Fhc = null_real !initialization: jauch
integer :: JX = null_int !initialization: jauch
integer :: IY = null_int !initialization: jauch
integer, dimension(:,:), pointer :: IZ => null()
integer, dimension(:,:), pointer :: JZ => null()
integer, dimension(:,:), pointer :: IU => null()
integer, dimension(:,:), pointer :: JU => null()
integer, dimension(:,:), pointer :: IV => null()
integer, dimension(:,:), pointer :: JV => null()
integer, dimension(:,:), pointer :: ICross => null()
integer, dimension(:,:), pointer :: JCross => null()
integer, dimension(:,:), pointer :: ILinkZ => null() !Joao Sobrinho
integer, dimension(:,:), pointer :: JLinkZ => null()
integer, dimension(:,:), pointer :: ILinkU => null()
integer, dimension(:,:), pointer :: JLinkU => null()
integer, dimension(:,:), pointer :: ILinkV => null()
integer, dimension(:,:), pointer :: JLinkV => null()
type (T_Size2D) :: MPI_Window
type (T_FatherGrid), pointer :: Next => null()
type (T_FatherGrid), pointer :: Prev => null()
end type T_FatherGrid
type T_Border
!Grid boundary
type(T_Polygon), pointer :: Polygon_ => null()
integer :: Type_ = null_int
end type T_Border
private :: T_Coef2D
type T_Coef2D
logical :: AllocateON = .false.
real(8), pointer, dimension(:) :: VECG => null()
real(8), pointer, dimension(:) :: VECW => null()
real, pointer, dimension(:,:) :: Results2D => null()
real, pointer, dimension(:,:) :: D => null()
real(8), pointer, dimension(:,:) :: E => null()
real , pointer, dimension(:,:) :: F => null()
real , pointer, dimension(:,:) :: Ti => null()
end type T_Coef2D
private :: T_Coef3D
type T_Coef3D
logical :: AllocateON = .false.
real(8), pointer, dimension(:) :: VECG => null()
real(8), pointer, dimension(:) :: VECW => null()
real, pointer, dimension(:,:,:) :: Results3D => null()
real, pointer, dimension(:,:,:) :: D => null()
real(8), pointer, dimension(:,:,:) :: E => null()
real, pointer, dimension(:,:,:) :: F => null()
real, pointer, dimension(:,:,:) :: Ti => null()
real, pointer, dimension(:,:,:) :: C => null()
real, pointer, dimension(:,:,:) :: G => null()
integer :: KLB = FillValueInt
integer :: KUB = FillValueInt
end type T_Coef3D
private :: T_DDecomp
type T_DDecomp
logical :: ON = .false.
logical :: Master = .false.
logical :: MasterOrSlave = .false.
logical :: Auto = .false.
integer :: Master_MPI_ID = null_int
integer :: Nslaves = 0
integer, dimension(:), pointer :: Slaves_MPI_ID => null()
type (T_Size2D), dimension(:), pointer :: Slaves_Size => null()
type (T_Size2D), dimension(:), pointer :: Slaves_Inner => null()
type (T_Size2D), dimension(:), pointer :: Slaves_Mapping => null()
type (T_Size2D), dimension(:), pointer :: Slaves_HaloMap => null()
integer :: MPI_ID = null_int
type (T_Size2D) :: Global
type (T_Size2D) :: Mapping
type (T_Size2D) :: Inner
type (T_Size2D) :: HaloMap
integer :: NeighbourSouth = null_int
integer :: NeighbourWest = null_int
integer :: NeighbourEast = null_int
integer :: NeighbourNorth = null_int
integer :: Halo_Points = null_int
integer, pointer, dimension(:,:) :: Interfaces => null()
integer :: NInterfaces = null_int
logical, dimension(1:4) :: OpenBordersON = .true.
character(PathLength) :: FilesListName = "DecomposedFiles.dat"
logical :: FilesListOpen = .false.
integer :: FilesListID = null_int
character(PathLength) :: ModelPath = null_str
type (T_Coef2D) :: Coef2D
type (T_Coef3D) :: Coef3D
logical :: AutomaticLines = .false.
end type T_DDecomp
type T_HorizontalGrid
integer :: InstanceID = null_int !initialization: jauch - or should be set to 0 (zero)?
!Former Bathymetry
real, pointer, dimension(: ) :: XX => null()
real, pointer, dimension(: ) :: YY => null()
real :: Xorig = null_real
real :: Yorig = null_real
real :: Latitude = null_real
real :: Longitude = null_real
real :: GRID_ANGLE = null_real
integer :: ProjType = null_int
real :: SP1 = null_real
real :: SP2 = null_real
real :: Easting = null_real !initialization: jauch
real :: Northing = null_real !initialization: jauch
integer :: Datum = null_int !ellipsoid
logical :: UseLambert = .FALSE.
integer :: CoordType = null_int
integer :: ZoneLong = null_int
integer :: ZoneLat = null_int
integer, dimension(2) :: Grid_Zone
integer, dimension(:, :), allocatable :: IWD_connections_U
integer, dimension(:, :), allocatable :: IWD_connections_V
integer, dimension(:, :), allocatable :: IWD_connections_Z
real, pointer, dimension(:) :: IWD_Distances_U => null()
real, pointer, dimension(:) :: IWD_Distances_V => null()
real, pointer, dimension(:) :: IWD_Distances_Z => null()
logical :: UsedIWD_2Way = .false.
integer :: IWD_Nodes_Z = null_int
integer :: IWD_Nodes_U = null_int
integer :: IWD_Nodes_V = null_int
type(T_Compute) :: Compute
type (T_FatherGrid), pointer :: FirstFatherGrid => null()
type (T_FatherGrid), pointer :: LastFatherGrid => null()
!Grid boundary
type(T_Border), pointer :: GridBorderCart => null()
type(T_Border), pointer :: GridBorderCoord => null()
type(T_Border), pointer :: GridOutBorderCoord => null()
type(T_Border), pointer :: GridBorderAlongGrid => null()
type(T_Border), pointer :: GridOutBorderCart => null()
!Distances (DXX... DVY)
real, dimension(:, :), pointer :: DXX => null()
real, dimension(:, :), pointer :: DYY => null()
real, dimension(:, :), pointer :: DZX => null()
real, dimension(:, :), pointer :: DZY => null()
real, dimension(:, :), pointer :: DUX => null()
real, dimension(:, :), pointer :: DUY => null()
real, dimension(:, :), pointer :: DVX => null()
real, dimension(:, :), pointer :: DVY => null()
real, dimension(:, :), pointer :: XX_IE => null()
real, dimension(:, :), pointer :: YY_IE => null()
real, dimension(:, :), pointer :: XX_AlongGrid => null()
real, dimension(:, :), pointer :: YY_AlongGrid => null()
real, dimension(:, :), pointer :: LatitudeConn => null()
real, dimension(:, :), pointer :: LongitudeConn => null()
logical :: ReadCartCorners = .false.
real, dimension(:, :), pointer :: RotationX => null()
real, dimension(:, :), pointer :: RotationY => null()
logical :: CornersXYInput = .false.
logical :: Distortion = .false.
logical :: RegularRotation = .false.
integer, dimension(:,:), pointer :: DefineCellsMap => null()
integer, dimension(:,:), pointer :: DefineFacesUMap => null()
integer, dimension(:,:), pointer :: DefineFacesVMap => null()
integer, dimension(:,:), pointer :: DefineCrossMap => null()
logical :: NotDefinedCells = .false.
!Latitude, Longitude
real, dimension(:, :), pointer :: LatitudeZ => null()
real, dimension(:, :), pointer :: LongitudeZ => null()
!Coriolis Factor
real, dimension(:, :), pointer :: F => null()
!Area
real, dimension(:, :), pointer :: GridCellArea => null()
!1D aux
real, dimension(: ), pointer :: XX1D_Aux => null()
real, dimension(: ), pointer :: YY1D_Aux => null()
!Other
type (T_Size2D) :: Size
type (T_Size2D) :: WorkSize
type (T_Size2D) :: GlobalWorkSize
character(PathLength) :: FileName = null_str
type(T_DDecomp) :: DDecomp
type (T_BorderLimits) :: BorderLimits
!Instances
integer :: ObjHDF5 = 0
integer :: ObjEnterData = 0
integer :: ObjEnterData2 = 0
type(T_Polygon), pointer :: AuxPolygon
type (T_HorizontalGrid), pointer :: GlobalGrid
type (T_HorizontalGrid), pointer :: Next => null()
end type T_HorizontalGrid
!Global Module Variables
type (T_HorizontalGrid), pointer :: FirstHorizontalGrid => null()
type (T_HorizontalGrid), pointer :: Me => null()
contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR CONSTRUCTOR
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine ConstructHorizontalGridV1(HorizontalGridID, DataFile, &
MPI_ID, MasterID, LastSlaveID, ModelPath, STAT)
!Arguments-------------------------------------------------------------
integer :: HorizontalGridID
character(len=*), optional :: DataFile
integer, optional, intent(IN) :: MPI_ID
integer, optional, intent(IN) :: MasterID
integer, optional, intent(IN) :: LastSlaveID
character(len=*), optional, intent(IN) :: ModelPath
integer, optional, intent(OUT) :: STAT
!Local-----------------------------------------------------------------
integer :: STAT_, ready_
integer :: STAT_CALL
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mHorizontalGrid_)) then
nullify (FirstHorizontalGrid)
call RegisterModule (mHorizontalGrid_)
endif
call Ready(HorizontalGridID, ready_)
cd2 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
nullify (Me%FirstFatherGrid)
nullify (Me%LastFatherGrid )
!Reads data file
if (.not. present(DataFile)) then
call ReadFileName('IN_BATIM', Me%FileName, Message = "Grid File", STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructHorizontalGridV1 - HorizontalGrid - ERR01'
else
Me%FileName = DataFile
endif
if (present(MPI_ID)) then
Me%DDecomp%MPI_ID = MPI_ID
endif
if (present(MasterID)) then
Me%DDecomp%Master_MPI_ID = MasterID
if (MasterID == MPI_ID) then
Me%DDecomp%Master = .true.
endif
endif
if (present(LastSlaveID)) then
if (MPI_ID > null_int .and. MasterID > null_int) then
if (LastSlaveID > null_int) then
Me%DDecomp%ON = .true.
Me%DDecomp%MasterOrSlave = .true.
Me%DDecomp%Nslaves = LastSlaveID - MasterID
if (present(ModelPath)) then
Me%DDecomp%ModelPath = ModelPath
else
Me%DDecomp%ModelPath = null_str
endif
endif
endif
endif
!Construct the variable common to all module
call ConstructGlobalVariables
call GenerateGrid
!Returns ID
HorizontalGridID = Me%InstanceID
STAT_ = SUCCESS_
else
stop 'HorizontalGrid - ConstructHorizontalGridV1 - ERR99'
end if cd2
if (present(STAT)) STAT = STAT_
end subroutine ConstructHorizontalGridV1
!--------------------------------------------------------------------------
subroutine ConstructHorizontalGridV2(HorizontalGridID, LatitudeConn, LongitudeConn, &
XX, YY, Xorig, Yorig, Latitude, Longitude, &
ILB, IUB, JLB, JUB, STAT)
!Arguments-------------------------------------------------------------
integer :: HorizontalGridID
real, dimension(:, :), pointer, optional :: LatitudeConn, LongitudeConn
real, dimension(: ), pointer :: XX, YY
real , optional :: Xorig, Yorig
real :: Latitude, Longitude
integer :: ILB, IUB, JLB, JUB
integer, optional, intent(OUT) :: STAT
!Local-----------------------------------------------------------------
real, dimension(:, :), pointer :: LatitudeConn_, LongitudeConn_
real :: Xorig_, Yorig_
integer :: STAT_, ready_
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mHorizontalGrid_)) then
nullify (FirstHorizontalGrid)
call RegisterModule (mHorizontalGrid_)
endif
call Ready(HorizontalGridID, ready_)
cd2 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
nullify (Me%FirstFatherGrid)
nullify (Me%LastFatherGrid )
if (present(LongitudeConn)) then
LongitudeConn_ => LongitudeConn
else
nullify(LongitudeConn_)
endif
if (present(LatitudeConn)) then
LatitudeConn_ => LatitudeConn
else
nullify(LatitudeConn_)
endif
if (present(Xorig)) then
Xorig_ = Xorig
else
Xorig_ = 0.
endif
if (present(Yorig)) then
Yorig_ = Yorig
else
Yorig_ = 0.
endif
call ConstructGlobalVariablesV1(LatitudeConn_, LongitudeConn_, XX, YY, &
Xorig_, Yorig_, Latitude, Longitude, &
ILB, IUB, JLB, JUB)
call GenerateGrid
!Returns ID
HorizontalGridID = Me%InstanceID
STAT_ = SUCCESS_
else
stop 'HorizontalGrid - ConstructHorizontalGridV2 - ERR99'
end if cd2
if (present(STAT)) STAT = STAT_
end subroutine ConstructHorizontalGridV2
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine ConstructHorizontalGridV3(HorizontalGridID, HDF5ID, STAT)
!Arguments-------------------------------------------------------------
integer :: HorizontalGridID
integer :: HDF5ID
integer, optional, intent(OUT) :: STAT
!Local-----------------------------------------------------------------
integer :: STAT_, ready_, nUsers
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
!Assures nullification of the global variable
if (.not. ModuleIsRegistered(mHorizontalGrid_)) then
nullify (FirstHorizontalGrid)
call RegisterModule (mHorizontalGrid_)
endif
call Ready(HorizontalGridID, ready_)
cd2 : if (ready_ .EQ. OFF_ERR_) then
call AllocateInstance
!Associates Horizontal Grid
Me%ObjHDF5 = AssociateInstance (mHDF5_, HDF5ID)
nullify (Me%FirstFatherGrid)
nullify (Me%LastFatherGrid )
call ConstructGlobalVariablesV2()
!Associates Horizontal Grid
nUsers = DeassociateInstance (mHDF5_, HDF5ID)
if (nUsers == 0) stop 'HorizontalGrid - ConstructHorizontalGridV3 - ERR10'
call GenerateGrid
!Returns ID
HorizontalGridID = Me%InstanceID
STAT_ = SUCCESS_
else
stop 'HorizontalGrid - ConstructHorizontalGridV3 - ERR20'
end if cd2
if (present(STAT)) STAT = STAT_
end subroutine ConstructHorizontalGridV3
!--------------------------------------------------------------------------
subroutine GenerateGrid
!Constructs XX_IE, YY_IE, LatitudeZ and LongitudeZ
call Mercator
if (.not. Me%CornersXYInput) then
!Constructs XX_Z, YY_Z, XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross
call GridPointsLocation1D
endif
!Constructs XX2D_Z, YY2D_Z, XX2D_U, YY2D_U, XX2D_V, YY2D_V, XX2D_Cross, YY2D_Cross
call GridPointsLocation2D
!Check grid border type
call CheckGridBorder
!Computes DXX, DYY, DZX, DZY, DUX, DUY, DVX, DVY, XX_AlongGrid, YY_AlongGrid
call ComputeDistances
!Computes RotationX, RotationY
call ComputeRotation
!Computes Area
call ComputeGridCellArea
!Computes Coriolis
call ComputeCoriolis
!Defines the grid border polygon
call DefineBorderPolygons
!Intialization of domain decomposition procedure
!call ConstructDDecomp
end subroutine GenerateGrid
subroutine ConstructDDecomp
#ifdef _USE_MPI
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
character(len = PathLength) :: DDFile
character(len = StringLength) :: Message
integer :: STAT_CALL, i, iflag
integer :: AuxHalo
!type (T_Size2D) :: Mapping
logical :: Exist
integer, pointer, dimension(:) :: Aux1D
integer :: Source, Destination
integer :: iSize
integer, save :: Precision
integer :: status(MPI_STATUS_SIZE)
!----------------------------------------------------------------------
Me%DDecomp%Global%JLB = Me%GlobalWorkSize%JLB
Me%DDecomp%Global%JUB = Me%GlobalWorkSize%JUB
Me%DDecomp%Global%ILB = Me%GlobalWorkSize%ILB
Me%DDecomp%Global%IUB = Me%GlobalWorkSize%IUB
allocate(Me%DDecomp%Slaves_MPI_ID(Me%DDecomp%Nslaves))
do i=1, Me%DDecomp%Nslaves
Me%DDecomp%Slaves_MPI_ID(i) = Me%DDecomp%Master_MPI_ID + i
enddo
! ---> ASCII file used to define domain decomposition
Message ='ASCII file used to define domain decomposition.'
Message = trim(Message)
!call ReadFileName('D_DECOMP', DDFile, Message = Message, STAT = STAT_CALL)
call GetData(DDFile,Me%ObjEnterData, iflag, &
keyword = 'D_DECOMP', &
default = null_str, &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, INTERNAL_, "ConstructDDecomp - Hydrodynamic - ERR10")
ii: if (iflag == 0) then
Me%DDecomp%Auto = .true.
else ii
inquire(FILE = DDFile, EXIST = Exist)
iE: if (Exist) then
!open()
call ConstructEnterData(Me%ObjEnterData2, DDFile, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, INTERNAL_, "ConstructDDecomp - Hydrodynamic - ERR20")
write(*,*) "Read from file Domain Decomposition mapping"
write(*,*) "Present MPI ID =", Me%DDecomp%MPI_ID
write(*,*) "Master MPI_ID = ", Me%DDecomp%Master_MPI_ID
write(*,*) "Number of domain Slaves = ",Me%DDecomp%Nslaves
do i=1, Me%DDecomp%Nslaves
write(*,*) "ID of Slave number ",i, "is =", Me%DDecomp%Slaves_MPI_ID(i)
enddo
call OptionsDDecomp
call KillEnterData(Me%ObjEnterData2, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) &
call SetError(FATAL_, INTERNAL_, "ConstructDDecomp - Hydrodynamic - ERR30")
else iE
write(*,*) 'keyword D_DECOMP do not exit file =',trim(DDFile)
write(*,*) 'Automatic decomposition will be assumed'
Me%DDecomp%Auto = .true.
endif iE
endif ii
if (Me%DDecomp%MasterOrSlave) then
if (Me%DDecomp%Auto) then
call AutomaticDDecomp
endif
if (Me%DDecomp%NeighbourSouth /= null_int) then
AuxHalo = Me%DDecomp%Halo_Points
!1 - South; 2 - North; 3 - West; 4 - East
Me%DDecomp%OpenBordersON(1) = .false.
else
AuxHalo = 0.
endif
Me%DDecomp%HaloMap%ILB = Me%DDecomp%Mapping%ILB - AuxHalo
Me%WorkSize%ILB = 1
Me%DDecomp%Inner%ILB = Me%WorkSize%ILB + AuxHalo
if (Me%DDecomp%NeighbourNorth /= null_int) then
AuxHalo = Me%DDecomp%Halo_Points
!1 - South; 2 - North; 3 - West; 4 - East
Me%DDecomp%OpenBordersON(2) = .false.
else
AuxHalo = 0.
endif
Me%DDecomp%HaloMap%IUB = Me%DDecomp%Mapping%IUB + AuxHalo
Me%WorkSize%IUB = Me%DDecomp%HaloMap%IUB - Me%DDecomp%HaloMap%ILB + 1
Me%DDecomp%Inner%IUB = Me%WorkSize%IUB - AuxHalo
if (Me%DDecomp%HaloMap%IUB - Me%DDecomp%HaloMap%ILB /= &
Me%WorkSize%IUB - Me%WorkSize%ILB) then
write(*,*) "Decomposition domains is inconsistent with the grid data input - Lines "
write(*,*) 'WorkSize ILB,IUB, =',Me%WorkSize%ILB,Me%WorkSize%IUB
write(*,*) 'HaloMap ILB,IUB, =',Me%DDecomp%HaloMap%ILB,Me%DDecomp%HaloMap%IUB
stop "ConstructDDecomp - ModuleHorizontalGrid - ERR40"
endif
if (Me%DDecomp%NeighbourWest /= null_int) then
AuxHalo = Me%DDecomp%Halo_Points
!1 - South; 2 - North; 3 - West; 4 - East
Me%DDecomp%OpenBordersON(3) = .false.
else
AuxHalo = 0.
endif
Me%DDecomp%HaloMap%JLB = Me%DDecomp%Mapping%JLB - AuxHalo
Me%WorkSize%JLB = 1
Me%DDecomp%Inner%JLB = Me%WorkSize%JLB + AuxHalo
if (Me%DDecomp%NeighbourEast /= null_int) then
AuxHalo = Me%DDecomp%Halo_Points
!1 - South; 2 - North; 3 - West; 4 - East
Me%DDecomp%OpenBordersON(4) = .false.
else
AuxHalo = 0.
endif
Me%DDecomp%HaloMap%JUB = Me%DDecomp%Mapping%JUB + AuxHalo
Me%WorkSize%JUB = Me%DDecomp%HaloMap%JUB - Me%DDecomp%HaloMap%JLB + 1
Me%DDecomp%Inner%JUB = Me%WorkSize%JUB - AuxHalo
if (Me%DDecomp%HaloMap%JUB - Me%DDecomp%HaloMap%JLB /= &
Me%WorkSize%JUB - Me%WorkSize%JLB) then
write(*,*) "Decomposition domains is inconsistent with the grid data input - Columns "
write(*,*) 'WorkSize JLB,JUB, =',Me%WorkSize%JLB,Me%WorkSize%JUB
write(*,*) 'HaloMap JLB,JUB, =',Me%DDecomp%HaloMap%JLB,Me%DDecomp%HaloMap%JUB
stop "ConstructDDecomp - ModuleHorizontalGrid - ERR50"
endif
if (Me%DDecomp%Master) then
!if (Me%DDecomp%Nslaves /= 7) write(*,*) 'e44'
!allocate(Me%DDecomp%Slaves_MPI_ID (Me%DDecomp%Nslaves))
allocate(Me%DDecomp%Slaves_Inner (Me%DDecomp%Nslaves))
allocate(Me%DDecomp%Slaves_Size (Me%DDecomp%Nslaves))
allocate(Me%DDecomp%Slaves_Mapping(Me%DDecomp%Nslaves))
allocate(Me%DDecomp%Slaves_HaloMap(Me%DDecomp%Nslaves))
endif
!PCL - Master slave mapping
iSize = 16
allocate(Aux1D(iSize))
if (Me%DDecomp%Master) then
do i=1, Me%DDecomp%Nslaves
Precision = MPI_INTEGER
Source = Me%DDecomp%Slaves_MPI_ID(i)
write(*,*) 'mpi_receive from to ', i, Me%DDecomp%MPI_ID
call MPI_Recv (Aux1D(1:iSize), iSize, Precision, Source, 30001, MPI_COMM_WORLD, status, STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop "ConstructDDecomp - ModuleHorizontalGrid - ERR60"
Me%DDecomp%Slaves_Inner (i)%ILB = Aux1D(1)
Me%DDecomp%Slaves_Inner (i)%IUB = Aux1D(2)
Me%DDecomp%Slaves_Inner (i)%JLB = Aux1D(3)
Me%DDecomp%Slaves_Inner (i)%JUB = Aux1D(4)
Me%DDecomp%Slaves_HaloMap(i)%ILB = Aux1D(5)
Me%DDecomp%Slaves_HaloMap(i)%IUB = Aux1D(6)
Me%DDecomp%Slaves_HaloMap(i)%JLB = Aux1D(7)
Me%DDecomp%Slaves_HaloMap(i)%JUB = Aux1D(8)
Me%DDecomp%Slaves_Size (i)%ILB = Aux1D(9)
Me%DDecomp%Slaves_Size (i)%IUB = Aux1D(10)
Me%DDecomp%Slaves_Size (i)%JLB = Aux1D(11)
Me%DDecomp%Slaves_Size (i)%JUB = Aux1D(12)
Me%DDecomp%Slaves_Mapping(i)%ILB = Aux1D(13)
Me%DDecomp%Slaves_Mapping(i)%IUB = Aux1D(14)
Me%DDecomp%Slaves_Mapping(i)%JLB = Aux1D(15)
Me%DDecomp%Slaves_Mapping(i)%JUB = Aux1D(16)
enddo
else
Aux1D(1) = Me%DDecomp%Inner%ILB
Aux1D(2) = Me%DDecomp%Inner%IUB
Aux1D(3) = Me%DDecomp%Inner%JLB
Aux1D(4) = Me%DDecomp%Inner%JUB
Aux1D(5) = Me%DDecomp%HaloMap%ILB
Aux1D(6) = Me%DDecomp%HaloMap%IUB
Aux1D(7) = Me%DDecomp%HaloMap%JLB
Aux1D(8) = Me%DDecomp%HaloMap%JUB
Aux1D(9) = Me%WorkSize%ILB
Aux1D(10)= Me%WorkSize%IUB
Aux1D(11)= Me%WorkSize%JLB
Aux1D(12)= Me%WorkSize%JUB
Aux1D(13) = Me%DDecomp%Mapping%ILB
Aux1D(14) = Me%DDecomp%Mapping%IUB
Aux1D(15) = Me%DDecomp%Mapping%JLB
Aux1D(16) = Me%DDecomp%Mapping%JUB
Precision = MPI_INTEGER
Destination = Me%DDecomp%Master_MPI_ID
write(*,*) 'mpi_send', Me%DDecomp%MPI_ID
call MPI_Send (Aux1D(1:iSize), iSize, Precision, Destination, 30001, MPI_COMM_WORLD, STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop "ConstructDDecomp - ModuleHorizontalGrid - ERR80"
endif
deallocate(Aux1d)
endif
#endif _USE_MPI
end subroutine ConstructDDecomp
!End----------------------------------------------------------------
#if _USE_MPI
subroutine OptionsDDecomp()
!Arguments------------------------------------------------------------
!Local----------------------------------------------------------------
character(len = StringLength), parameter :: BeginBlock1 ="<BeginSubDD>"
character(len = StringLength), parameter :: EndBlock1 ="<EndSubDD>"
character(len = StringLength), parameter :: BeginBlock2 ="<BeginInterfaceSN>"
character(len = StringLength), parameter :: EndBlock2 ="<EndInterfaceSN>"
character(len = StringLength), parameter :: BeginBlock3 ="<BeginInterfaceWE>"
character(len = StringLength), parameter :: EndBlock3 ="<EndInterfaceWE>"
integer, dimension(:), allocatable :: Aux1D
logical :: BlockFound
integer :: ClientNumber, STAT_CALL, iflag
integer :: in, line, FirstLine, LastLine, i, ii, jj
integer :: SN_N_Interfaces, WE_N_Interfaces, MPI_ID
logical :: MissMatchID
!Begin----------------------------------------------------------------
call RewindBuffer(Me%ObjEnterData2, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR120'
Me%DDecomp%Auto = .false.
iSl: do i =1, Me%DDecomp%Nslaves + 1
!Searches sub-domains blocks
call ExtractBlockFromBuffer (Me%ObjEnterData2, ClientNumber, &
BeginBlock1, EndBlock1, &
BlockFound, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR130'
if (.not. BlockFound ) then
Me%DDecomp%Auto = .true.
write(*,*) 'OptionsDDecomp - ModuleHorizontalGrid - WRN140'
write(*,*) 'Domain Decomposition in automatic mode'
exit
endif
call GetData(Value = MPI_ID, &
EnterDataID = Me%ObjEnterData2, &
flag = iflag, &
keyword = 'MPI_ID', &
SearchType = FromBlock, &
ClientModule = 'ModuleHorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR150'
if (iflag == 0 ) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR160'
MissMatchID = .true.
do ii =1, Me%DDecomp%Nslaves
if (Me%DDecomp%Slaves_MPI_ID(ii) == MPI_ID) then
MissMatchID = .false.
endif
enddo
if (Me%DDecomp%Master_MPI_ID == MPI_ID) then
MissMatchID = .false.
endif
if (MissMatchID) then
write(*,*) 'MPI ID of present Domain', Me%DDecomp%MPI_ID
write(*,*) 'Domain -', MPI_ID, ' is not one of decomposition domains'
write(*,*) "All MPI_ID - Slaves"
do ii =1, Me%DDecomp%Nslaves
write(*,*) "MPI_ID =", Me%DDecomp%Slaves_MPI_ID(ii)
enddo
write(*,*) "MPI_ID Master=",Me%DDecomp%Master_MPI_ID
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR165'
endif
if (MPI_ID == Me%DDecomp%MPI_ID) then
Me%DDecomp%MasterOrSlave = .true.
allocate(Aux1D(1:2))
call GetData(Vector = Aux1D, &
EnterDataID = Me%ObjEnterData2, &
flag = iflag, &
keyword = 'ILB_IUB', &
SearchType = FromBlock, &
ClientModule = 'ModuleHorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_ .and. STAT_CALL /= KEYWORD_NOT_FOUND_ERR_) then
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR170'
endif
if (iflag /= 2) then
if (iflag == 0) then
Aux1D(1) = Me%DDecomp%Global%ILB
Aux1D(2) = Me%DDecomp%Global%IUB
else
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR180'
endif
endif
Me%DDecomp%Mapping%ILB = Aux1D(1)
Me%DDecomp%Mapping%IUB = Aux1D(2)
call GetData(Vector = Aux1D, &
EnterDataID = Me%ObjEnterData2, &
flag = iflag, &
keyword = 'JLB_JUB', &
SearchType = FromBlock, &
ClientModule = 'ModuleHorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_ .and. STAT_CALL /= KEYWORD_NOT_FOUND_ERR_) then
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR190'
endif
if (iflag /= 2) then
if (iflag == 0) then
Aux1D(1) = Me%DDecomp%Global%JLB
Aux1D(2) = Me%DDecomp%Global%JUB
else
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR200'
endif
endif
Me%DDecomp%Mapping%JLB = Aux1D(1)
Me%DDecomp%Mapping%JUB = Aux1D(2)
deallocate(Aux1D)
exit
endif
enddo iSl
call Block_Unlock(Me%ObjEnterData2, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR210'
iAuto: if (Me%DDecomp%Auto) then
call GetData(Value = Me%DDecomp%AutomaticLines, &
EnterDataID = Me%ObjEnterData2, &
flag = iflag, &
keyword = 'AUTOMATIC_LINES', &
SearchType = FromFile, &
default = .false., &
ClientModule = 'ModuleHorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR230'
else
call GetData(Value = Me%DDecomp%NInterfaces, &
EnterDataID = Me%ObjEnterData2, &
flag = iflag, &
keyword = 'INTERFACES_NUMBER', &
SearchType = FromFile, &
ClientModule = 'ModuleHorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR230'
allocate(Me%DDecomp%Interfaces(Me%DDecomp%NInterfaces,3))
allocate(Aux1D(1:2))
call RewindBuffer(Me%ObjEnterData2, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR220'
!Searches sub-domains blocks
call ExtractBlockFromBuffer (Me%ObjEnterData2, ClientNumber, &
BeginBlock2, EndBlock2, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR240'
if (.not. BlockFound) then
SN_N_Interfaces = 0
endif
in = 0
if (BlockFound) then
SN_N_Interfaces = LastLine - FirstLine - 1
if (LastLine == FirstLine) then
SN_N_Interfaces = 0
endif
if (SN_N_Interfaces > Me%DDecomp%NInterfaces) then
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR250'
endif
do line = FirstLine + 1, LastLine - 1
call GetData(Vector = Aux1D, &
EnterDataID = Me%ObjEnterData2, &
flag = iflag, &
SearchType = FromBlock, &
ClientModule = 'ModuleHorizontalGrid', &
Buffer_Line = line, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR260'
if (iflag /= 2 ) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR270'
in = in + 1
do jj = 1, 2
MissMatchID = .true.
do ii =1, Me%DDecomp%Nslaves
if (Me%DDecomp%Slaves_MPI_ID(ii) == Aux1D(jj)) then
MissMatchID = .false.
endif
enddo
if (Me%DDecomp%Master_MPI_ID == Aux1D(jj)) then
MissMatchID = .false.
endif
if (MissMatchID) then
write(*,*) 'Domain -', Aux1D(jj), ' is not one of decomposition domains'
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR275'
endif
enddo
Me%DDecomp%Interfaces(in,1) = Aux1D(1)
Me%DDecomp%Interfaces(in,2) = Aux1D(2)
Me%DDecomp%Interfaces(in,3) = SouthNorth_
if (Me%DDecomp%MPI_ID == Aux1D(1)) then
Me%DDecomp%NeighbourNorth = Aux1D(2)
endif
if (Me%DDecomp%MPI_ID == Aux1D(2)) then
Me%DDecomp%NeighbourSouth = Aux1D(1)
endif
enddo
endif
call RewindBuffer(Me%ObjEnterData2, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR275'
!Searches sub-domains blocks
call ExtractBlockFromBuffer (Me%ObjEnterData2, ClientNumber, &
BeginBlock3, EndBlock3, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR280'
if (.not. BlockFound) then
WE_N_Interfaces = 0
endif
if (BlockFound) then
WE_N_Interfaces = LastLine - FirstLine - 1
if (LastLine == FirstLine) then
WE_N_Interfaces = 0
endif
if (SN_N_Interfaces + WE_N_Interfaces /= Me%DDecomp%NInterfaces) then
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR290'
endif
do line = FirstLine + 1, LastLine - 1
call GetData(Vector = Aux1D, &
EnterDataID = Me%ObjEnterData2, &
flag = iflag, &
SearchType = FromBlock, &
ClientModule = 'ModuleHorizontalGrid', &
Buffer_Line = line, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR300'
if (iflag /= 2 ) stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR310'
in = in + 1
do jj = 1, 2
MissMatchID = .true.
do ii =1, Me%DDecomp%Nslaves
if (Me%DDecomp%Slaves_MPI_ID(ii) == Aux1D(jj)) then
MissMatchID = .false.
endif
enddo
if (Me%DDecomp%Master_MPI_ID == Aux1D(jj)) then
MissMatchID = .false.
endif
if (MissMatchID) then
write(*,*) 'Domain -', Aux1D(jj), ' is not one of decomposition domains'
stop 'OptionsDDecomp - ModuleHorizontalGrid - ERR315'
endif
enddo
Me%DDecomp%Interfaces(in,1) = Aux1D(1)
Me%DDecomp%Interfaces(in,2) = Aux1D(2)
Me%DDecomp%Interfaces(in,3) = WestEast_
if (Me%DDecomp%MPI_ID == Aux1D(1)) then
Me%DDecomp%NeighbourEast = Aux1D(2)
endif
if (Me%DDecomp%MPI_ID == Aux1D(2)) then
Me%DDecomp%NeighbourWest = Aux1D(1)
endif
enddo
endif
deallocate(Aux1D)
endif iAuto
end subroutine OptionsDDecomp
!End------------------------------------------------------------------
subroutine AutomaticDDecomp()
!Arguments------------------------------------------------------------
!Local----------------------------------------------------------------
!Begin----------------------------------------------------------------
write(*,*) 'halo_points', Me%DDecomp%Halo_Points
if (Me%DDecomp%Global%IUB > Me%DDecomp%Global%JUB .or. Me%DDecomp%AutomaticLines) then
call AutomaticDDecompLines ()
else
call AutomaticDDecompColumns()
endif
end subroutine AutomaticDDecomp
!End------------------------------------------------------------------
subroutine AutomaticDDecompLines()
!Arguments------------------------------------------------------------
!Local----------------------------------------------------------------
integer :: i, iSouth, iNorth
integer :: iub_map, ilb_map, ILB, IUB, ND, is
real :: IDD
!Begin----------------------------------------------------------------
!In automatic mode MOHID slices the domains along the matrix lines (ILB-IUB)
Me%DDecomp%Mapping%JLB = Me%DDecomp%Global%JLB
Me%DDecomp%Mapping%JUB = Me%DDecomp%Global%JUB
ILB = Me%DDecomp%Global%ILB
IUB = Me%DDecomp%Global%IUB
ND = Me%DDecomp%Nslaves + 1
IDD = real (IUB - ILB + 1) / real(ND)
is = Me%DDecomp%Master_MPI_ID
do i=1, Me%DDecomp%Nslaves + 1
if (is + i - 1 == Me%DDecomp%MPI_ID) then
ilb_map = ILB + int(real(i-1)*IDD)
iub_map = ILB - 1 + int(real(i )*IDD)
exit
endif
enddo
Me%DDecomp%Mapping%ILB = ilb_map
Me%DDecomp%Mapping%IUB = iub_map
write(*,*) 'Limits ', Me%DDecomp%MPI_ID, iub_map, ilb_map
Me%DDecomp%NInterfaces = Me%DDecomp%Nslaves
allocate(Me%DDecomp%Interfaces(Me%DDecomp%NInterfaces,3))
do i = 1, Me%DDecomp%NInterfaces
iSouth = Me%DDecomp%Master_MPI_ID + i - 1
iNorth = Me%DDecomp%Master_MPI_ID + i
Me%DDecomp%Interfaces(i,1) = iSouth
Me%DDecomp%Interfaces(i,2) = iNorth
Me%DDecomp%Interfaces(i,3) = SouthNorth_
write(*,*) 'Interface ', i, Me%DDecomp%Interfaces(i,1), Me%DDecomp%Interfaces(i,2)
if (Me%DDecomp%MPI_ID == iSouth) then
Me%DDecomp%NeighbourNorth = iNorth
endif
if (Me%DDecomp%MPI_ID == iNorth) then
Me%DDecomp%NeighbourSouth = iSouth
endif
enddo
end subroutine AutomaticDDecompLines
!End------------------------------------------------------------------
subroutine AutomaticDDecompColumns()
!Arguments------------------------------------------------------------
!Local----------------------------------------------------------------
integer :: i, jWest, jEast
integer :: jub_map, jlb_map, JLB, JUB, ND, is
real :: IDD
!Begin----------------------------------------------------------------
!In automatic mode MOHID slices the domains along the matrix columns (JLB-JUB)
Me%DDecomp%Mapping%ILB = Me%DDecomp%Global%ILB
Me%DDecomp%Mapping%IUB = Me%DDecomp%Global%IUB
JLB = Me%DDecomp%Global%JLB
JUB = Me%DDecomp%Global%JUB
ND = Me%DDecomp%Nslaves + 1
IDD = real (JUB - JLB + 1) / real(ND)
is = Me%DDecomp%Master_MPI_ID
do i=1, Me%DDecomp%Nslaves + 1
if (is + i - 1 == Me%DDecomp%MPI_ID) then
jlb_map = JLB + int(real(i-1)*IDD)
jub_map = JLB - 1 + int(real(i )*IDD)
exit
endif
enddo
Me%DDecomp%Mapping%JLB = jlb_map
Me%DDecomp%Mapping%JUB = jub_map
write(*,*) 'Limits ', Me%DDecomp%MPI_ID, jub_map, jlb_map
Me%DDecomp%NInterfaces = Me%DDecomp%Nslaves
allocate(Me%DDecomp%Interfaces(Me%DDecomp%NInterfaces,3))
do i = 1, Me%DDecomp%NInterfaces
jWest = Me%DDecomp%Master_MPI_ID + i - 1
jEast = Me%DDecomp%Master_MPI_ID + i
Me%DDecomp%Interfaces(i,1) = jWest
Me%DDecomp%Interfaces(i,2) = jEast
Me%DDecomp%Interfaces(i,3) = WestEast_
write(*,*) 'Interface ', i, Me%DDecomp%Interfaces(i,1), Me%DDecomp%Interfaces(i,2)
if (Me%DDecomp%MPI_ID == jWest) then
Me%DDecomp%NeighbourEast = jEast
endif
if (Me%DDecomp%MPI_ID == jEast) then
Me%DDecomp%NeighbourWest = jWest
endif
enddo
end subroutine AutomaticDDecompColumns
!End------------------------------------
#endif _USE_MPI
subroutine AllocateInstance
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
type (T_HorizontalGrid), pointer :: NewObjHorizontalGrid
type (T_HorizontalGrid), pointer :: PreviousObjHorizontalGrid
!Allocates new instance
allocate (NewObjHorizontalGrid)
nullify (NewObjHorizontalGrid%Next)
!Insert New Instance into list and makes Current point to it
if (.not. associated(FirstHorizontalGrid)) then
FirstHorizontalGrid => NewObjHorizontalGrid
Me => NewObjHorizontalGrid
else
PreviousObjHorizontalGrid => FirstHorizontalGrid
Me => FirstHorizontalGrid%Next
do while (associated(Me))
PreviousObjHorizontalGrid => Me
Me => Me%Next
enddo
Me => NewObjHorizontalGrid
PreviousObjHorizontalGrid%Next => NewObjHorizontalGrid
endif
Me%InstanceID = RegisterNewInstance (mHORIZONTALGRID_)
end subroutine AllocateInstance
!--------------------------------------------------------------------------
subroutine ConstructFatherGridLocation(HorizontalGridID, HorizontalGridFatherID, &
GridID, OkCross, OkZ, OkU, OkV, Window, STAT)
!Arguments-------------------------------------------------------------
integer :: HorizontalGridID
integer :: HorizontalGridFatherID
integer, optional, intent (IN) :: GridID
logical, optional, intent (IN) :: OkCross, OkZ, OkU, OkV
type (T_Size2D), optional :: Window
integer, optional, intent (OUT) :: STAT
!Local-----------------------------------------------------------------
type (T_HorizontalGrid), pointer :: ObjHorizontalGridFather
type (T_FatherGrid), pointer :: NewFatherGrid
integer :: STAT_, ready_, GridID_
logical :: OkZ_, OkU_, OkV_, OkCross_
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
call Ready(HorizontalGridID, ready_)
cd1 : if (ready_ .EQ. IDLE_ERR_) then
if (present(GridID)) then
GridID_ = GridID
else
GridID_ = LargeScaleModel_
endif
if (present(OkCross)) then
OkCross_ = OkCross
else
OkCross_ = .false.
endif
if (present(OkZ)) then
OkZ_ = OkZ
else
OkZ_ = .true.
endif
if (present(OkU)) then
OkU_ = OkU
else
OkU_ = .true.
endif
if (present(OkV)) then
OkV_ = OkV
else
OkV_ = .true.
endif
allocate(NewFatherGrid)
nullify(ObjHorizontalGridFather)
call LocateObjFather (ObjHorizontalGridFather, HorizontalGridFatherID)
if (present(Window)) then
NewFatherGrid%MPI_Window = Window
else
NewFatherGrid%MPI_Window = ObjHorizontalGridFather%WorkSize
endif
if (ObjHorizontalGridFather%CornersXYInput) then
call ConstructNewFatherGrid2D (ObjHorizontalGridFather, NewFatherGrid, GridID_, &
OkZ_, OkU_, OkV_, OkCross_)
else
call ConstructNewFatherGrid1D (ObjHorizontalGridFather, NewFatherGrid, GridID_, &
OkZ_, OkU_, OkV_, OkCross_)
endif
if (Me%CoordType == SIMPLE_GEOG_ .and. .not. Me%ReadCartCorners .and. Me%ProjType == PAULO_PROJECTION_) then
if (ObjHorizontalGridFather%Longitude /= Me%Longitude) then
write(*,*) 'LONGITUDE (in the bathymetry file) must be the same in Father and Son Domains'
stop 'ConstructFatherGridLocation - ModuleHorizontalGrid - ERR10'
endif
if (ObjHorizontalGridFather%Latitude /= Me%Latitude) then
write(*,*) 'LATITUDE (in the bathymetry file) must be the same in Father and Son Domains'
stop 'ConstructFatherGridLocation - ModuleHorizontalGrid - ERR20'
endif
endif
call Add_FatherGrid (NewFatherGrid)
STAT_ = SUCCESS_
else
STAT_ = ready_
end if cd1
if (present(STAT)) STAT = STAT_
end subroutine ConstructFatherGridLocation
!--------------------------------------------------------------------------
subroutine GetSonWindow(HorizontalGridID, HorizontalGridFatherID, Window, STAT)
!Arguments-------------------------------------------------------------
integer :: HorizontalGridID
integer :: HorizontalGridFatherID
type (T_Size2D), intent (OUT) :: Window
integer, optional, intent (OUT) :: STAT
!Local-----------------------------------------------------------------
type (T_HorizontalGrid), pointer :: ObjHorizontalGridFather
integer :: STAT_, ready_
real, dimension(:,:), pointer :: XX2D, YY2D
real, dimension(:,:), pointer :: WindowInXY
integer, dimension(:,:), pointer :: WindowOutJI
logical :: WindowWithData
real :: West, East, South, North
integer :: ILB, IUB, JLB, JUB
!------------------------------------------------------------------------
STAT_ = UNKNOWN_
call Ready(HorizontalGridID, ready_)
cd1 : if (ready_ .EQ. IDLE_ERR_) then
nullify(ObjHorizontalGridFather)
call LocateObjFather (ObjHorizontalGridFather, HorizontalGridFatherID)
allocate (WindowInXY (1:2,1:2))
allocate (WindowOutJI(1:2,1:2))
WindowOutJI(:,:) = null_int
WindowInXY (:,:) = null_real
West = Me%GridOutBorderCoord%Polygon_%Limits%Left
East = Me%GridOutBorderCoord%Polygon_%Limits%Right
South = Me%GridOutBorderCoord%Polygon_%Limits%Bottom
North = Me%GridOutBorderCoord%Polygon_%Limits%Top
WindowInXY(2,1) = South
WindowInXY(2,2) = North
WindowInXY(1,1) = West
WindowInXY(1,2) = East
ILB = ObjHorizontalGridFather%WorkSize%ILB
IUB = ObjHorizontalGridFather%WorkSize%IUB
JLB = ObjHorizontalGridFather%WorkSize%JLB
JUB = ObjHorizontalGridFather%WorkSize%JUB
if (Me%CoordType == SIMPLE_GEOG_ .or. Me%CoordType == GEOG_) then
XX2D => ObjHorizontalGridFather%LongitudeConn
YY2D => ObjHorizontalGridFather%LatitudeConn
else
XX2D => ObjHorizontalGridFather%XX_IE
YY2D => ObjHorizontalGridFather%YY_IE
endif
call ArrayPolygonWindow(XX = XX2D, &
YY = YY2D, &
WIn = WindowInXY, &
ILB = ILB, &
IUB = IUB+1, &
JLB = JLB, &
JUB = JUB+1, &
WOut = WindowOutJI, &
WindowWithData = WindowWithData)
if (.not.WindowWithData) then
write(*,*) 'Father grid do not intersect Nested model'
stop 'GetSonWindow - ModuleHoriuzontalGrid - ERR20'
endif
Window%ILB = max(WindowOutJI(1,1) - 1,ILB)
Window%IUB = min(WindowOutJI(1,2) + 1,IUB)
Window%JLB = max(WindowOutJI(2,1) - 1,JLB)
Window%JUB = min(WindowOutJI(2,2) + 1,JUB)
deallocate (WindowInXY )
deallocate (WindowOutJI)
STAT_ = SUCCESS_
else
STAT_ = ready_
end if cd1
if (present(STAT)) STAT = STAT_
end subroutine GetSonWindow
!--------------------------------------------------------------------------
subroutine GetTwoWayAux (HorizontalGridID, IWD_Connections_U, IWD_Connections_V, IWD_Connections_Z, &
IWD_Distances_U, IWD_Distances_V, IWD_Distances_Z, IWD_Nodes_Z, IWD_Nodes_U, &
IWD_Nodes_V, STAT)
!Arguments-------------------------------------------------------------
integer :: HorizontalGridID, IWD_Nodes_Z, IWD_Nodes_U, IWD_Nodes_V
integer, optional, intent (OUT) :: STAT
integer, dimension(:,:), pointer :: IWD_Connections_U, IWD_Connections_V, IWD_Connections_Z
real, dimension(: ), pointer :: IWD_Distances_U, IWD_Distances_V, IWD_Distances_Z
!Local-----------------------------------------------------------------
integer :: STAT_, ready_
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
call Ready(HorizontalGridID, ready_)
if ((ready_ .EQ. IDLE_ERR_ ) .OR. &
(ready_ .EQ. READ_LOCK_ERR_)) then
IWD_Connections_U => Me%IWD_Connections_U
call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
IWD_Connections_V => Me%IWD_Connections_V
call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
IWD_Connections_Z => Me%IWD_Connections_Z
call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
IWD_Distances_U => Me%IWD_Distances_U
call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
IWD_Distances_V => Me%IWD_Distances_V
call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
IWD_Distances_Z => Me%IWD_Distances_Z
call Read_Lock(mHORIZONTALGRID_, Me%InstanceID)
IWD_Nodes_U = Me%IWD_Nodes_U
IWD_Nodes_V = Me%IWD_Nodes_V
IWD_Nodes_Z = Me%IWD_Nodes_Z
STAT_ = SUCCESS_
else
STAT_ = ready_
end if
if (present(STAT)) &
STAT = STAT_
end subroutine GetTwoWayAux
!--------------------------------------------------------------------------
subroutine UnGetTwoWayAux(HorizontalGridID, IWD_Connections_U, IWD_Connections_V, IWD_Connections_Z, &
IWD_Distances_U, IWD_Distances_V, IWD_Distances_Z, STAT)
!Arguments-------------------------------------------------------------
integer :: HorizontalGridID
integer, optional, intent (OUT) :: STAT
integer, dimension(:,:), pointer :: IWD_Connections_U, IWD_Connections_V, IWD_Connections_Z
real, dimension(: ), pointer :: IWD_Distances_U, IWD_Distances_V, IWD_Distances_Z
!Local-----------------------------------------------------------------
integer :: STAT_, ready_
!----------------------------------------------------------------------
STAT_ = UNKNOWN_
call Ready(HorizontalGridID, ready_)
if ((ready_ .EQ. IDLE_ERR_ ) .OR. &
(ready_ .EQ. READ_LOCK_ERR_)) then
nullify(IWD_Connections_U)
call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UnGetTwoWayAux - ERR01")
nullify(IWD_Connections_V)
call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UnGetTwoWayAux - ERR02")
nullify(IWD_Connections_Z)
call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UnGetTwoWayAux - ERR03")
nullify(IWD_Distances_U)
call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UnGetTwoWayAux - ERR04")
nullify(IWD_Distances_V)
call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UnGetTwoWayAux - ERR05")
nullify(IWD_Distances_Z)
call Read_UnLock(mHORIZONTALGRID_, Me%InstanceID, "UnGetTwoWayAux - ERR06")
STAT_ = SUCCESS_
else
STAT_ = ready_
end if
if (present(STAT)) &
STAT = STAT_
end subroutine UnGetTwoWayAux
!--------------------------------------------------------------------------
subroutine ConstructNewFatherGrid1D(ObjHorizontalGridFather, NewFatherGrid, GridID, &
OkZ, OkU, OkV, OkCross, CheckRotation)
!Arguments-------------------------------------------------------------
type (T_HorizontalGrid), pointer :: ObjHorizontalGridFather
type (T_FatherGrid) , pointer :: NewFatherGrid
integer, intent (IN) :: GridID
logical, intent (IN) :: OkZ, OkU, OkV, OkCross
logical, optional, intent (IN) :: CheckRotation
!Local-----------------------------------------------------------------
real, dimension(:,:), pointer :: XX_Z, YY_Z, XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross
integer, dimension(:,:), pointer :: IZ, JZ, IU, JU, IV, JV, ICross, JCross
real :: Xorig , Yorig , Rotation
integer :: ILB, IUB, JLB, JUB, i, j
integer :: ILBwork, IUBwork, JLBwork, JUBwork
logical :: CheckRotation_
integer :: ILBFAther, IUBFather, JLBFather, JUBFather, U, V
!----------------------------------------------------------------------
!if (Me%CornersXYInput) stop 'ConstructNewFatherGrid1D - ModuleHoriuzontalGrid - ERR01'
if (present(CheckRotation)) then
CheckRotation_ = CheckRotation
else
CheckRotation_ = .true.
endif
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
ILBwork = Me%WorkSize%ILB
IUBwork = Me%WorkSize%IUB
JLBwork = Me%WorkSize%JLB
JUBwork = Me%WorkSize%JUB
NewFatherGrid%GridID = GridID
NewFatherGrid%OkZ = OkZ
NewFatherGrid%OkU = OkU
NewFatherGrid%OkV = OkV
NewFatherGrid%OkCross = OkCross
NewFatherGrid%CornersXYInput = .false.
nullify (NewFatherGrid%XX_Z )
nullify (NewFatherGrid%YY_Z )
nullify (NewFatherGrid%XX_U )
nullify (NewFatherGrid%YY_U )
nullify (NewFatherGrid%XX_V )
nullify (NewFatherGrid%YY_V )
nullify (NewFatherGrid%Next )
nullify (NewFatherGrid%Prev )
!Gets the bounds from the Father
ILBFAther = NewFatherGrid%MPI_Window%ILB
IUBFather = NewFatherGrid%MPI_Window%IUB
JLBFather = NewFatherGrid%MPI_Window%JLB
JUBFather = NewFatherGrid%MPI_Window%JUB
NewFatherGrid%MPI_Window%ILB = -null_int
NewFatherGrid%MPI_Window%IUB = null_int
NewFatherGrid%MPI_Window%JLB = -null_int
NewFatherGrid%MPI_Window%JUB = null_int
!Compute points location ib the father grid
if (NewFatherGrid%OkZ) then
!Allocates Variables
allocate (NewFatherGrid%XX_Z (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%YY_Z (ILB:IUB, JLB:JUB))
!Initialize Values
NewFatherGrid%XX_Z (:,:) = FillValueReal
NewFatherGrid%YY_Z (:,:) = FillValueReal
allocate (NewFatherGrid%IZ (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JZ (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%ILinkZ (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JLinkZ (ILB:IUB, JLB:JUB))
NewFatherGrid%IZ (:,:) = FillValueInt
NewFatherGrid%JZ (:,:) = FillValueInt
NewFatherGrid%ILinkZ (:,:) = FillValueInt
NewFatherGrid%JLinkZ (:,:) = FillValueInt
XX_Z => NewFatherGrid%XX_Z
YY_Z => NewFatherGrid%YY_Z
IZ => NewFatherGrid%IZ
JZ => NewFatherGrid%JZ
endif
if (NewFatherGrid%OkU) then
allocate (NewFatherGrid%XX_U (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%YY_U (ILB:IUB, JLB:JUB))
NewFatherGrid%XX_U (:,:) = FillValueReal
NewFatherGrid%YY_U (:,:) = FillValueReal
allocate (NewFatherGrid%IU (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JU (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%ILinkU (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JLinkU (ILB:IUB, JLB:JUB))
NewFatherGrid%IU (:,:) = FillValueInt
NewFatherGrid%JU (:,:) = FillValueInt
NewFatherGrid%ILinkU (:,:) = FillValueInt
NewFatherGrid%JLinkU (:,:) = FillValueInt
XX_U => NewFatherGrid%XX_U
YY_U => NewFatherGrid%YY_U
IU => NewFatherGrid%IU
JU => NewFatherGrid%JU
endif
if (NewFatherGrid%OkV) then
allocate (NewFatherGrid%XX_V (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%YY_V (ILB:IUB, JLB:JUB))
NewFatherGrid%XX_V (:,:) = FillValueReal
NewFatherGrid%YY_V (:,:) = FillValueReal
allocate (NewFatherGrid%IV (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JV (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%ILinkV (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JLinkV (ILB:IUB, JLB:JUB))
NewFatherGrid%IV (:,:) = FillValueInt
NewFatherGrid%JV (:,:) = FillValueInt
NewFatherGrid%ILinkV (:,:) = FillValueInt
NewFatherGrid%JLinkV (:,:) = FillValueInt
YY_V => NewFatherGrid%YY_V
XX_V => NewFatherGrid%XX_V
IV => NewFatherGrid%IV
JV => NewFatherGrid%JV
endif
if (NewFatherGrid%OkCross) then
allocate (NewFatherGrid%XX_Cross (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%YY_Cross (ILB:IUB, JLB:JUB))
NewFatherGrid%XX_Cross (:,:) = FillValueReal
NewFatherGrid%YY_Cross (:,:) = FillValueReal
allocate (NewFatherGrid%ICross (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JCross (ILB:IUB, JLB:JUB))
NewFatherGrid%ICross (:,:) = FillValueInt
NewFatherGrid%JCross (:,:) = FillValueInt
YY_Cross => NewFatherGrid%YY_Cross
XX_Cross => NewFatherGrid%XX_Cross
ICross => NewFatherGrid%ICross
JCross => NewFatherGrid%JCross
endif
if (Me%CornersXYInput) then
Xorig = - ObjHorizontalGridFather%Xorig
Yorig = - ObjHorizontalGridFather%Yorig
Rotation = - ObjHorizontalGridFather%Grid_Angle
else
Xorig = Me%Xorig - ObjHorizontalGridFather%Xorig
Yorig = Me%Yorig - ObjHorizontalGridFather%Yorig
Rotation = Me%Grid_Angle - ObjHorizontalGridFather%Grid_Angle
if (Rotation /= 0..and. CheckRotation_) then
write(*,*) 'The submodels do not work with referentials with different rotations'
write(*,*) 'ConstructNewFatherGrid1D - HorizontalGrid - WARN10'
endif
call RODAXY(0., 0., -Me%Grid_Angle, Xorig, Yorig)
if (Xorig < 0 .and. Yorig < 0) then
write(*,*) 'The submodels must be locate inside the father domain'
write(*,*) 'ConstructNewFatherGrid1D - HorizontalGrid - WARN20'
endif
endif
!Compute points location ib the father grid
if (NewFatherGrid%OkZ) then
!Rotation of the points
do3 : do j = JLBwork, JUBwork
do4 : do i = ILBwork, IUBwork
if (Me%CornersXYInput) then
if (Me%DefineCellsMap(i,j) == 0) then
IZ(i, j) = FillValueInt
JZ(i, j) = FillValueInt
cycle
endif
XX_Z(i, j) = Me%Compute%XX2D_Z(i, j)
YY_Z(i, j) = Me%Compute%YY2D_Z(i, j)
else
XX_Z(i, j) = Me%Compute%XX_Z(j)
YY_Z(i, j) = Me%Compute%YY_Z(i)
endif
call RODAXY(Xorig, Yorig, Rotation, XX_Z(i, j), YY_Z(i, j))
call LocateCell (ObjHorizontalGridFather%Compute%XX_Z, &
ObjHorizontalGridFather%Compute%YY_Z, &
XX_Z(i, j), YY_Z(i, j), &
ILBFAther, IUBFather, JLBFather, JUBFather, &
IZ(i, j), JZ(i, j))
call LocateCell_2(ObjHorizontalGridFather%Compute%XX_U, &
ObjHorizontalGridFather%Compute%YY_V, &
XX_Z(i, j), YY_Z(i, j), &
ILBFAther, IUBFather, JLBFather, JUBFather, &
NewFatherGrid%ILinkZ(i, j), NewFatherGrid%JLinkZ(i, j))
!Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IZ(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IZ(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JZ(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JZ(i, j))
end do do4
end do do3
nullify(XX_Z , YY_Z )
nullify(IZ , JZ )
endif
if (NewFatherGrid%OkU) then
!Rotation of the points
do1 : do j = JLBwork, JUBwork+1
do2 : do i = ILBwork, IUBwork
if (Me%CornersXYInput) then
if (Me%DefineFacesUMap(i,j) == 0) then
IU(i, j) = FillValueInt
JU(i, j) = FillValueInt
cycle
endif
XX_U(i, j) = Me%Compute%XX2D_U(i, j)
YY_U(i, j) = Me%Compute%YY2D_U(i, j)
else
XX_U(i, j) = Me%Compute%XX_U(j)
YY_U(i, j) = Me%Compute%YY_U(i)
endif
call RODAXY(Xorig, Yorig, Rotation, XX_U(i, j), YY_U(i, j))
call LocateCell (ObjHorizontalGridFather%Compute%XX_U, &
ObjHorizontalGridFather%Compute%YY_U, &
XX_U(i, j), YY_U(i, j), &
ILBFAther, IUBFather, JLBFather, JUBFather + 1, &
IU(i, j), JU(i, j))
call LocateCell_2(ObjHorizontalGridFather%Compute%XX_Z, &
ObjHorizontalGridFather%Compute%YY_V, &
XX_U(i, j), YY_U(i, j), &
ILBFAther, IUBFather, JLBFather, JUBFather, &
NewFatherGrid%ILinkU(i, j), NewFatherGrid%JLinkU(i, j), U = U)
!Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IU(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IU(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JU(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JU(i, j))
end do do2
end do do1
nullify(XX_U , YY_U )
nullify(IU , JU )
endif
if (NewFatherGrid%OkV) then
!Rotation of the points
do5: do j = JLBwork, JUBwork
do6: do i = ILBwork, IUBwork+1
if (Me%CornersXYInput) then
if (Me%DefineFacesVMap(i,j) == 0) then
IV(i, j) = FillValueInt
JV(i, j) = FillValueInt
cycle
endif
XX_V(i, j) = Me%Compute%XX2D_V(i, j)
YY_V(i, j) = Me%Compute%YY2D_V(i, j)
else
XX_V(i, j) = Me%Compute%XX_V(j)
YY_V(i, j) = Me%Compute%YY_V(i)
endif
call RODAXY(Xorig, Yorig, Rotation, XX_V(i, j), YY_V(i, j))
call LocateCell (ObjHorizontalGridFather%Compute%XX_V, &
ObjHorizontalGridFather%Compute%YY_V, &
XX_V(i, j), YY_V(i, j), &
ILBFAther, IUBFather + 1, JLBFather, JUBFather, &
IV(i, j), JV(i, j))
call LocateCell_2(ObjHorizontalGridFather%Compute%XX_U, &
ObjHorizontalGridFather%Compute%YY_Z, &
XX_V(i, j), YY_V(i, j), &
ILBFAther, IUBFather, JLBFather, JUBFather, &
NewFatherGrid%ILinkV(i, j), NewFatherGrid%JLinkV(i, j), V = V)
!MPI_Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IV(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IV(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JV(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JV(i, j))
end do do6
end do do5
nullify(XX_V , YY_V )
nullify(IV , JV )
endif
if (NewFatherGrid%OkCross) then
!Rotation of the points
do7: do j = JLBwork, JUBwork
do8: do i = ILBwork, IUBwork
if (Me%CornersXYInput) then
if (Me%DefineCrossMap(i,j) == 0) then
ICross(i, j) = FillValueInt
JCross(i, j) = FillValueInt
cycle
endif
XX_Cross(i, j) = Me%XX_IE(i, j)
YY_Cross(i, j) = Me%YY_IE(i, j)
else
XX_Cross(i, j) = Me%Compute%XX_Z(j)
YY_Cross(i, j) = Me%Compute%YY_Z(i)
endif
call RODAXY(Xorig, Yorig, Rotation, XX_Cross(i, j), YY_Cross(i, j))
call LocateCell (ObjHorizontalGridFather%Compute%XX_Cross, &
ObjHorizontalGridFather%Compute%YY_Cross, &
XX_Cross(i, j), YY_Cross(i, j), &
ILBFAther + 1, IUBFather + 1, JLBFather, JUBFather, &
ICross(i, j), JCross(i, j))
!MPI_Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, ICross(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, ICross(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JCross(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JCross(i, j))
end do do8
end do do7
nullify(XX_Cross, YY_Cross)
nullify(ICross , JCross )
endif
!Increase MPI Window Size by one (Fluxes at the upper boundary)
NewFatherGrid%MPI_Window%IUB = NewFatherGrid%MPI_Window%IUB + 1
NewFatherGrid%MPI_Window%JUB = NewFatherGrid%MPI_Window%JUB + 1
end subroutine ConstructNewFatherGrid1D
!--------------------------------------------------------------------------
subroutine ConstructNewFatherGrid2D(ObjHorizontalGridFather, NewFatherGrid, GridID, &
OkZ, OkU, OkV, OkCross)
!Arguments-------------------------------------------------------------
type (T_HorizontalGrid), pointer :: ObjHorizontalGridFather
type (T_FatherGrid) , pointer :: NewFatherGrid
integer, intent (IN) :: GridID
logical, intent (IN) :: OkZ, OkU, OkV, OkCross
!Local-----------------------------------------------------------------
real :: XPoint, YPoint
integer, dimension(:,:), pointer :: IZ, JZ, IU, JU, IV, JV, &
ICross, JCross, DefinedPoint
integer :: ILB, IUB, JLB, JUB, i, j
integer :: ILBwork, IUBwork, JLBwork, JUBwork
integer :: ILBFAther, IUBFather, JLBFather, JUBFather
!----------------------------------------------------------------------
if (.not. Me%CornersXYInput) stop 'ConstructNewFatherGrid2D - ModuleHorizontalGrid - ERR01'
ILB = Me%Size%ILB
IUB = Me%Size%IUB
JLB = Me%Size%JLB
JUB = Me%Size%JUB
ILBwork = Me%WorkSize%ILB
IUBwork = Me%WorkSize%IUB
JLBwork = Me%WorkSize%JLB
JUBwork = Me%WorkSize%JUB
NewFatherGrid%GridID = GridID
NewFatherGrid%OkZ = OkZ
NewFatherGrid%OkU = OkU
NewFatherGrid%OkV = OkV
NewFatherGrid%OkCross = OkCross
!Gets the bounds from the Father
ILBFAther = NewFatherGrid%MPI_Window%ILB
IUBFather = NewFatherGrid%MPI_Window%IUB
JLBFather = NewFatherGrid%MPI_Window%JLB
JUBFather = NewFatherGrid%MPI_Window%JUB
NewFatherGrid%MPI_Window%ILB = -null_int
NewFatherGrid%MPI_Window%IUB = null_int
NewFatherGrid%MPI_Window%JLB = -null_int
NewFatherGrid%MPI_Window%JUB = null_int
NewFatherGrid%CornersXYInput = .true.
nullify (NewFatherGrid%Next )
nullify (NewFatherGrid%Prev )
write(*,*) 'The submodels in curvilinear grids must have identical distortion'
write(*,*) 'of the coarser grid'
write(*,*) 'ConstructNewFatherGrid2D - HorizontalGrid - WARN01'
!Compute points location ib the father grid
if (NewFatherGrid%OkZ) then
!Allocates Variables
allocate (NewFatherGrid%IZ (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JZ (ILB:IUB, JLB:JUB))
!Initialize Values
NewFatherGrid%IZ (:,:) = FillValueInt
NewFatherGrid%JZ (:,:) = FillValueInt
IZ => NewFatherGrid%IZ
JZ => NewFatherGrid%JZ
NewFatherGrid%XX_Z => Me%Compute%XX2D_Z
NewFatherGrid%YY_Z => Me%Compute%YY2D_Z
DefinedPoint => ObjHorizontalGridFather%DefineCellsMap
endif
if (NewFatherGrid%OkU) then
allocate (NewFatherGrid%IU (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JU (ILB:IUB, JLB:JUB))
NewFatherGrid%IU (:,:) = FillValueInt
NewFatherGrid%JU (:,:) = FillValueInt
IU => NewFatherGrid%IU
JU => NewFatherGrid%JU
NewFatherGrid%XX_U => Me%Compute%XX2D_U
NewFatherGrid%YY_U => Me%Compute%YY2D_U
DefinedPoint => ObjHorizontalGridFather%DefineFacesUMap
endif
if (NewFatherGrid%OkV) then
allocate (NewFatherGrid%IV (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JV (ILB:IUB, JLB:JUB))
NewFatherGrid%IV (:,:) = FillValueInt
NewFatherGrid%JV (:,:) = FillValueInt
IV => NewFatherGrid%IV
JV => NewFatherGrid%JV
NewFatherGrid%XX_V => Me%Compute%XX2D_V
NewFatherGrid%YY_V => Me%Compute%YY2D_V
DefinedPoint => ObjHorizontalGridFather%DefineFacesVMap
endif
if (NewFatherGrid%OkCross) then
allocate (NewFatherGrid%ICross (ILB:IUB, JLB:JUB))
allocate (NewFatherGrid%JCross (ILB:IUB, JLB:JUB))
NewFatherGrid%ICross (:,:) = FillValueInt
NewFatherGrid%JCross (:,:) = FillValueInt
ICross => NewFatherGrid%ICross
JCross => NewFatherGrid%JCross
NewFatherGrid%XX_Cross => Me%XX_IE
NewFatherGrid%YY_Cross => Me%YY_IE
DefinedPoint => ObjHorizontalGridFather%DefineCrossMap
endif
!Compute points location ib the father grid
if (NewFatherGrid%OkZ) then
!Rotation of the points
do3 : do j = JLBwork, JUBwork
do4 : do i = ILBwork, IUBwork
if (Me%CornersXYInput) then
XPoint = Me%Compute%XX2D_Z(i, j)
YPoint = Me%Compute%YY2D_Z(i, j)
else
XPoint = Me%Compute%XX_Z(j)
YPoint = Me%Compute%YY_Z(i)
endif
call LocateCellPolygons(ObjHorizontalGridFather%Compute%XX2D_Z, &
ObjHorizontalGridFather%Compute%YY2D_Z, &
XPoint, YPoint, DefinedPoint, &
ILBFAther, IUBFather, JLBFather, JUBFather, &
IZ(i, j), JZ(i, j))
!Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IZ(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IZ(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JZ(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JZ(i, j))
end do do4
end do do3
nullify(IZ , JZ )
endif
if (NewFatherGrid%OkU) then
!Rotation of the points
do1 : do j = JLBwork, JUBwork+1
do2 : do i = ILBwork, IUBwork
call LocateCellPolygons(ObjHorizontalGridFather%Compute%XX2D_U, &
ObjHorizontalGridFather%Compute%YY2D_U, &
Me%Compute%XX2D_U(i, j), &
Me%Compute%YY2D_U(i, j), &
DefinedPoint, &
ILBFAther, IUBFather, JLBFather, JUBFather + 1, &
IU(i, j), JU(i, j))
!Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IU(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IU(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JU(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JU(i, j))
end do do2
end do do1
nullify(IU , JU )
endif
if (NewFatherGrid%OkV) then
!Rotation of the points
do5: do j = JLBwork, JUBwork
do6: do i = ILBwork, IUBwork+1
call LocateCellPolygons(ObjHorizontalGridFather%Compute%XX2D_V, &
ObjHorizontalGridFather%Compute%YY2D_V, &
Me%Compute%XX2D_V(i, j), &
Me%Compute%YY2D_V(i, j), &
DefinedPoint, &
ILBFAther, IUBFather + 1, JLBFather, JUBFather, &
IV(i, j), JV(i, j))
!MPI_Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, IV(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, IV(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JV(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JV(i, j))
end do do6
end do do5
nullify(IV , JV )
endif
if (NewFatherGrid%OkCross) then
!Rotation of the points
do7: do j = JLBwork, JUBwork
do8: do i = ILBwork, IUBwork
call LocateCellPolygons(ObjHorizontalGridFather%XX_IE, &
ObjHorizontalGridFather%YY_IE, &
Me%XX_IE(i, j), &
Me%YY_IE(i, j), &
DefinedPoint, &
ILBFAther + 1, IUBFather + 1, JLBFather, JUBFather,&
ICross(i, j), JCross(i, j))
!MPI_Window
NewFatherGrid%MPI_Window%ILB = min(NewFatherGrid%MPI_Window%ILB, ICross(i, j))
NewFatherGrid%MPI_Window%IUB = max(NewFatherGrid%MPI_Window%IUB, ICross(i, j))
NewFatherGrid%MPI_Window%JLB = min(NewFatherGrid%MPI_Window%JLB, JCross(i, j))
NewFatherGrid%MPI_Window%JUB = max(NewFatherGrid%MPI_Window%JUB, JCross(i, j))
end do do8
end do do7
nullify(ICross , JCross )
endif
!Increase MPI Window Size by one (Fluxes at the upper boundary)
NewFatherGrid%MPI_Window%IUB = NewFatherGrid%MPI_Window%IUB + 1
NewFatherGrid%MPI_Window%JUB = NewFatherGrid%MPI_Window%JUB + 1
end subroutine ConstructNewFatherGrid2D
!-------------------------------------------------------------------------------------
!>@author Joao Sobrinho Maretec
!>@Brief
!>Builds connection matrix between son and father grid on a IWD interpolation
!>@param[in] FatherID, SonID
subroutine ConstructIWDTwoWay(FatherID, SonID)
!Arguments--------------------------------------------------------------
type (T_HorizontalGrid), pointer :: ObjHorizontalFather
integer :: FatherID, SonID
!Local------------------------------------------------------------------
real :: FatherCenterX, FatherCenterY, DistanceToFather, SonCenterX, SonCenterY
real :: SearchRadious, MaxRatio
integer :: index, Nbr_Connections, minJ, minI, maxJ, maxI
integer :: i, j, i2, j2, ready_
integer, dimension (:, :), pointer :: FatherLinkI, FatherLinkJ
!-------------------------------------------------------------------------
index = 1
call Ready (SonID, ready_)
call LocateObjFather (ObjHorizontalFather, FatherID)
call DetermineMaxRatio(ObjHorizontalFather, MaxRatio)
FatherLinkI => Me%LastFatherGrid%ILinkZ
FatherLinkJ => Me%LastFatherGrid%JLinkZ
minJ = min(FatherLinkJ(1,1), FatherLinkJ(1, Me%Size%JUB - 1), &
FatherLinkJ(Me%Size%IUB - 1, 1), FatherLinkJ(Me%Size%IUB - 1, Me%Size%JUB - 1))
minI = min(FatherLinkI(1,1), FatherLinkI(1, Me%Size%JUB - 1), &
FatherLinkI(Me%Size%IUB - 1, 1), FatherLinkI(Me%Size%IUB - 1, Me%Size%JUB - 1))
maxJ = max(FatherLinkJ(1,1), FatherLinkJ(1, Me%Size%JUB - 1), &
FatherLinkJ(Me%Size%IUB - 1, 1), FatherLinkJ(Me%Size%IUB - 1, Me%Size%JUB - 1))
maxI = max(FatherLinkI(1,1), FatherLinkI(1, Me%Size%JUB - 1), &
FatherLinkI(Me%Size%IUB - 1, 1), FatherLinkI(Me%Size%IUB - 1, Me%Size%JUB - 1))
nullify (FatherLinkI)
nullify (FatherLinkJ)
!Uses the maxRatio to avoid allocating too few indexes. 2nd term is to account for search radious
Nbr_Connections = (maxJ - minJ + 2) * (maxI - minI + 2) * (MaxRatio + 4 * (int(sqrt(MaxRatio))) + 4)
!Vectorials and scalars
allocate (Me%IWD_connections_Z (Nbr_Connections, 4))
allocate (Me%IWD_Distances_Z (Nbr_Connections))
Me%IWD_connections_Z(:, :) = FillValueInt
Me%IWD_Distances_Z(:) = FillValueReal
allocate (Me%IWD_connections_U(Nbr_Connections, 4))
allocate (Me%IWD_Distances_U (Nbr_Connections))
Me%IWD_connections_U(:, :) = FillValueInt
Me%IWD_Distances_U(:) = FillValueReal
allocate (Me%IWD_connections_V (Nbr_Connections, 4))
allocate (Me%IWD_Distances_V (Nbr_Connections))
Me%IWD_connections_V(:, :) = FillValueInt
Me%IWD_Distances_V(:) = FillValueReal
call ConstructIWDVel(ObjHorizontalFather, minI, minJ, maxI, maxJ, MaxRatio)
!i and j -> father cell
!i2 and j2 -> son cell
do j = minJ, maxJ
do i = minI, maxI
!Find Father's cell center coordinates
FatherCenterX = (( ObjHorizontalFather%XX_IE(i, j ) + ObjHorizontalFather%XX_IE(i+1, j ))/2. + &
( ObjHorizontalFather%XX_IE(i, j+1) + ObjHorizontalFather%XX_IE(i+1, j+1))/2.)/2.
FatherCenterY = (( ObjHorizontalFather%YY_IE(i, j ) + ObjHorizontalFather%YY_IE(i+1, j ))/2. + &
( ObjHorizontalFather%YY_IE(i, j+1) + ObjHorizontalFather%YY_IE(i+1, j+1))/2.)/2.
SearchRadious = (1.01+(1/(Sqrt(MaxRatio)))) * Sqrt((FatherCenterX - ObjHorizontalFather%XX(j))**2 + &
(FatherCenterY - ObjHorizontalFather%YY(i))**2)
!Find and build matrix of correspondent son cells
do j2 = 1, Me%Size%JUB - 1
do i2 = 1, Me%Size%IUB - 1
SonCenterX = (( Me%XX_IE(i2, j2 ) + Me%XX_IE(i2+1, j2 ))/2. + &
( Me%XX_IE(i2, j2+1) + Me%XX_IE(i2+1, j2+1))/2.)/2.
SonCenterY = (( Me%YY_IE(i2, j2 ) + Me%YY_IE(i2+1, j2 ))/2. + &
( Me%YY_IE(i2, j2+1) + Me%YY_IE(i2+1, j2+1))/2.)/2.
DistanceToFather = Sqrt((SonCenterX - FatherCenterX)**2.0 + &
(SonCenterY - FatherCenterY)**2.0)
if (DistanceToFather <= SearchRadious) then
Me%IWD_connections_Z(index, 1) = i
Me%IWD_connections_Z(index, 2) = j
Me%IWD_connections_Z(index, 3) = i2
Me%IWD_connections_Z(index, 4) = j2
if (DistanceToFather == 0)then
Me%IWD_Distances_Z(index) = 1.e-5
else
Me%IWD_Distances_Z(index) = DistanceToFather
endif
index = index + 1
endif
enddo
enddo
enddo
enddo
Me%IWD_Nodes_Z = index - 1
end subroutine ConstructIWDTwoWay
!---------------------------------------------------------------------------
!>@author Joao Sobrinho Maretec
!>@Brief
!>Builds connection matrix for velocities cells between son and father grid on a IWD interpolation
!>@param[in] ObjHorizontalFather, IWD_Distances, IWD_connections, minI, minJ, maxI, maxJ, MaxRatio
subroutine ConstructIWDVel(ObjHorizontalFather, minI, minJ, maxI, maxJ, MaxRatio)
!Arguments--------------------------------------------------------------
type (T_HorizontalGrid), pointer :: ObjHorizontalFather
integer :: minJ, minI, maxJ, maxI
real :: MaxRatio
!Local------------------------------------------------------------------
real :: FatherCenterX_U, FatherCenterX_V, FatherCenterY_U, FatherCenterY_V, &
FatherCenterX_Z, FatherCenterY_Z, SonCenterX_U, SonCenterX_V, &
SonCenterY_U, SonCenterY_V, DistanceToFather_U, DistanceToFather_V, &
SearchRadious_U, SearchRadious_V
integer :: index_U, index_V, i, j, i2, j2
!-------------------------------------------------------------------------
index_U = 1
index_V = 1
do j = minJ, maxJ
do i = minI, maxI
!Find Father cell center for U cell and V cell
FatherCenterX_U = (ObjHorizontalFather%XX_IE(i, j ))/2. + (ObjHorizontalFather%XX_IE(i+1, j ))/2.
FatherCenterY_U = (ObjHorizontalFather%YY_IE(i, j ))/2. + (ObjHorizontalFather%YY_IE(i+1, j ))/2.
FatherCenterX_V = (ObjHorizontalFather%XX_IE(i, j ))/2. + (ObjHorizontalFather%XX_IE(i , j+1))/2.
FatherCenterY_V = (ObjHorizontalFather%YY_IE(i, j ))/2. + (ObjHorizontalFather%YY_IE(i , j+1))/2.
FatherCenterX_Z = (( ObjHorizontalFather%XX_IE(i, j ) + ObjHorizontalFather%XX_IE(i+1, j ))/2. + &
( ObjHorizontalFather%XX_IE(i, j+1) + ObjHorizontalFather%XX_IE(i+1, j+1))/2.)/2.
FatherCenterY_Z = (( ObjHorizontalFather%YY_IE(i, j ) + ObjHorizontalFather%YY_IE(i+1, j ))/2. + &
( ObjHorizontalFather%YY_IE(i, j+1) + ObjHorizontalFather%YY_IE(i+1, j+1))/2.)/2.
SearchRadious_U = (1.01+(1/(Sqrt(MaxRatio)))) * &
Sqrt((FatherCenterX_U - FatherCenterX_Z)**2 + &
(FatherCenterY_U - ObjHorizontalFather%YY_IE(i, j))**2)
SearchRadious_V = (1.01+(1/(Sqrt(MaxRatio)))) * &
Sqrt((FatherCenterX_V - ObjHorizontalFather%XX_IE(i, j))**2 + &
(FatherCenterY_V - FatherCenterY_Z)**2)
! Find and build matrix of correspondent son cells
do j2 = 1, Me%Size%JUB - 1
do i2 = 1, Me%Size%IUB - 1
SonCenterX_U = (Me%XX_IE(i2, j2 ))/2. + (Me%XX_IE(i2+1, j2 ))/2.
SonCenterY_U = (Me%YY_IE(i2, j2 ))/2. + (Me%YY_IE(i2+1, j2 ))/2.
SonCenterX_V = (Me%XX_IE(i2, j2 ))/2. + (Me%XX_IE(i2 , j2+1))/2.
SonCenterY_V = (Me%YY_IE(i2, j2 ))/2. + (Me%YY_IE(i2 , j2+1))/2.
DistanceToFather_U = Sqrt((SonCenterX_U - FatherCenterX_U)**2.0 + &
(SonCenterY_U - FatherCenterY_U)**2.0)
DistanceToFather_V = Sqrt((SonCenterX_V - FatherCenterX_V)**2.0 + &
(SonCenterY_V - FatherCenterY_V)**2.0)
if (DistanceToFather_U <= SearchRadious_U) then
Me%IWD_connections_U(index_U, 1) = i
Me%IWD_connections_U(index_U, 2) = j
Me%IWD_connections_U(index_U, 3) = i2
Me%IWD_connections_U(index_U, 4) = j2
if (DistanceToFather_U == 0)then
!The 0.001 is the reference distance. The if also avoids /0 in module functions
Me%IWD_Distances_U(index_U) = 1.e-5
else
Me%IWD_Distances_U(index_U) = DistanceToFather_U
endif
index_U = index_U + 1
endif
if (DistanceToFather_V <= SearchRadious_V) then
Me%IWD_connections_V(index_V, 1) = i
Me%IWD_connections_V(index_V, 2) = j
Me%IWD_connections_V(index_V, 3) = i2
Me%IWD_connections_V(index_V, 4) = j2
if (DistanceToFather_V == 0)then
!The 0.001 is the reference distance. The if also avoids /0 in module functions
Me%IWD_Distances_V(index_V) = 1.e-5
else
Me%IWD_Distances_V(index_V) = DistanceToFather_V
endif
index_V = index_V + 1
endif
enddo
enddo
enddo
enddo
Me%IWD_Nodes_U = index_U - 1
Me%IWD_Nodes_V = index_V - 1
end subroutine ConstructIWDVel
!---------------------------------------------------------------------------
!>@author Joao Sobrinho Maretec
!>@Brief
!>Computes Maximum grid ratio between father and son domains
!>@param[in] ObjHorizontalGridFather, Ratio
subroutine DetermineMaxRatio(ObjHorizontalGridFather, MaxRatio)
!Argumments---------------------------------------------------------------------
type (T_HorizontalGrid), pointer :: ObjHorizontalGridFather
real, intent(OUT) :: MaxRatio
!local--------------------------------------------------------------------------
integer :: i, j
real :: aux
!-------------------------------------------------------------------------------
MaxRatio = 1.0
aux = 1.0
do j = 1, Me%Size%JUB - 1
do i = 1, Me%Size%IUB - 1
aux = ObjHorizontalGridFather%GridCellArea(Me%LastFatherGrid%IZ(i, j), Me%LastFatherGrid%JZ(i, j)) / &
Me%GridCellArea(i, j)
if (aux > MaxRatio) then
MaxRatio = aux
endif
enddo
enddo
end subroutine DetermineMaxRatio
!--------------------------------------------------------------------------
! This subroutine adds a new grid to the father grid List
subroutine Add_FatherGrid(NewFatherGrid)
!Arguments-------------------------------------------------------------
type(T_FatherGrid), pointer :: NewFatherGrid
!----------------------------------------------------------------------
! Add to the WaterFatherGrid List a new father grid
if (.not. associated(Me%FirstFatherGrid)) then
Me%FirstFatherGrid => NewFatherGrid
Me%LastFatherGrid => NewFatherGrid
else
NewFatherGrid%Prev => Me%LastFatherGrid
Me%LastFatherGrid%Next => NewFatherGrid
Me%LastFatherGrid => NewFatherGrid
end if
!----------------------------------------------------------------------
end subroutine Add_FatherGrid
!--------------------------------------------------------------------------
subroutine Search_FatherGrid(NewFatherGrid, GridID)
!Arguments-------------------------------------------------------------
type(T_FatherGrid), pointer :: NewFatherGrid
integer, intent (IN) :: GridID
!Local-----------------------------------------------------------------
NewFatherGrid => Me%FirstFatherGrid
do1 : do while (associated(NewFatherGrid))
cd1 : if (NewFatherGrid%GridID == GridID) then
exit do1
else
NewFatherGrid => NewFatherGrid%Next
end if cd1
end do do1
if (.not.Associated(NewFatherGrid)) &
call SetError(FATAL_, INTERNAL_, 'Search_FatherGrid - HorizontalGrid - ERR01')
!----------------------------------------------------------------------
end subroutine Search_FatherGrid
!--------------------------------------------------------------------------
subroutine ConstructGlobalVariables
!Arguments-------------------------------------------------------------
!Local-----------------------------------------------------------------
integer :: STAT_CALL, ZoneLong
integer, dimension(2) :: AuxInt
real, dimension(2) :: AuxReal
real :: DX, DY, Aux, XY_Aux
integer :: flag, flag1, flag2
integer :: ClientNumber
logical :: BlockFound, ConstantSpacingX, ConstantSpacingY
integer :: FirstLine, LastLine, line, i, j, ii, jj, iflag
!----------------------------------------------------------------------
!Nullify T_Distances
nullify (Me%XX )
nullify (Me%YY )
nullify (Me%DXX)
nullify (Me%DYY)
nullify (Me%DZX)
nullify (Me%DZY)
nullify (Me%DUX)
nullify (Me%DUY)
nullify (Me%DVX)
nullify (Me%DVY)
nullify (Me%XX_IE)
nullify (Me%YY_IE)
nullify (Me%Compute%XX_Z )
nullify (Me%Compute%YY_Z )
nullify (Me%Compute%XX_U )
nullify (Me%Compute%YY_U )
nullify (Me%Compute%XX_V )
nullify (Me%Compute%YY_V )
nullify (Me%Compute%XX_Cross)
nullify (Me%Compute%YY_Cross)
nullify (Me%XX_AlongGrid)
nullify (Me%YY_AlongGrid)
!Nullify Other
nullify (Me%LatitudeZ )
nullify (Me%LongitudeZ )
nullify (Me%F )
nullify (Me%GridCellArea )
nullify (Me%LatitudeConn )
nullify (Me%LongitudeConn)
!Opens Data File
call ConstructEnterData(Me%ObjEnterData, Me%FileName, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR10'
!Reads ILB_IUB
call GetData(AuxInt,Me%ObjEnterData, flag, &
keyword = 'ILB_IUB', &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR20'
if (flag /= 2 ) stop 'ConstructGlobalVariables - HorizontalGrid - ERR30'
Me%GlobalWorkSize%ILB = AuxInt(1)
Me%GlobalWorkSize%IUB = AuxInt(2)
!Reads JLB_JUB
call GetData(AuxInt,Me%ObjEnterData, flag, &
keyword = 'JLB_JUB', &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR40'
if (flag /= 2 ) stop 'ConstructGlobalVariables - HorizontalGrid - ERR50'
Me%GlobalWorkSize%JLB = AuxInt(1)
Me%GlobalWorkSize%JUB = AuxInt(2)
call GetData(Value = Me%DDecomp%Halo_Points, &
EnterDataID = Me%ObjEnterData, &
flag = iflag, &
keyword = 'HALOPOINTS', &
SearchType = FromFile, &
ClientModule = 'HorizontalGrid', &
default = 3, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR440'
if (Me%DDecomp%Halo_Points < 2) then
write(*,*) 'HALOPOINTS must equal to 2 or bigger'
stop 'ConstructGlobalVariables - HorizontalGrid - ERR445'
endif
!Intialization of domain decomposition procedure
if (Me%DDecomp%ON) then
call ConstructDDecomp
endif
if (Me%DDecomp%MasterOrSlave) then
Me%Size%ILB = Me%WorkSize%ILB-1
Me%Size%IUB = Me%WorkSize%IUB+1
if (Me%GlobalWorkSize%ILB /= Me%DDecomp%Global%ILB) then
write(*,*) "Me%GlobalWorkSize%ILB",Me%GlobalWorkSize%ILB
write(*,*) "Me%DDecomp%Global%ILB",Me%DDecomp%Global%ILB
stop 'ConstructGlobalVariables - HorizontalGrid - ERR32'
endif
if (Me%GlobalWorkSize%IUB /= Me%DDecomp%Global%IUB) then
stop 'ConstructGlobalVariables - HorizontalGrid - ERR34'
endif
else
Me%Size%ILB = Me%GlobalWorkSize%ILB-1
Me%Size%IUB = Me%GlobalWorkSize%IUB+1
Me%WorkSize%ILB = Me%GlobalWorkSize%ILB
Me%WorkSize%IUB = Me%GlobalWorkSize%IUB
endif
if (Me%DDecomp%MasterOrSlave) then
Me%Size%JLB = Me%WorkSize%JLB-1
Me%Size%JUB = Me%WorkSize%JUB+1
if (Me%GlobalWorkSize%JLB /= Me%DDecomp%Global%JLB) then
stop 'ConstructGlobalVariables - HorizontalGrid - ERR52'
endif
if (Me%GlobalWorkSize%JUB /= Me%DDecomp%Global%JUB) then
stop 'ConstructGlobalVariables - HorizontalGrid - ERR54'
endif
else
Me%Size%JLB = Me%GlobalWorkSize%JLB-1
Me%Size%JUB = Me%GlobalWorkSize%JUB+1
Me%WorkSize%JLB = Me%GlobalWorkSize%JLB
Me%WorkSize%JUB = Me%GlobalWorkSize%JUB
endif
!Reads COORD_TIP
call GetData(Me%CoordType, Me%ObjEnterData, flag, &
keyword = 'COORD_TIP', &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR60'
if (.NOT. ((Me%CoordType == GEOG_ ) .OR. &
(Me%CoordType == UTM_ ) .OR. &
(Me%CoordType == MIL_PORT_ ) .OR. &
(Me%CoordType == SIMPLE_GEOG_) .OR. &
(Me%CoordType == GRID_COORD_ ) .OR. &
(Me%CoordType == CIRCULAR_ ) .OR. &
(Me%CoordType == NLRD_ ))) then
call SetError (FATAL_, INTERNAL_, "Invalid Coordinates")
endif
!Reads Latitude (Reference latitude. If externally specificed is the center of domain)
call GetData(Me%Latitude, Me%ObjEnterData, flag, &
keyword = 'LATITUDE', &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR70'
!Reads Longitude (Reference longitude. If externally specificed is the center of domain)
call GetData(Me%Longitude, Me%ObjEnterData, flag, &
keyword = 'LONGITUDE', &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR80'
call GetData(Me%Datum, Me%ObjEnterData, flag, &
keyword = 'DATUM', &
default = WGS_84_DATUM, &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR90'
!Valor metrico a acrescentar as coordenadas na direccao x
call GetData(Me%Easting, Me%ObjEnterData, flag, &
keyword = 'EASTING', &
default = 0.0, &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR100'
!Valor metrico a acrescentar as coordenadas na direccao x
call GetData(Me%Northing, Me%ObjEnterData, flag, &
keyword = 'NORTHING', &
default = 0.0, &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR110'
!Projection Type to convert to cartesian coordinates
call GetData(Me%ProjType, Me%ObjEnterData, flag, &
keyword = 'PROJ_TYPE', &
default = PAULO_PROJECTION_, &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR120'
!Lambert Conformal Conic Projection
call GetData(Me%SP1, Me%ObjEnterData, flag1, &
keyword = 'SP1', &
default = 30., &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR120'
!Lambert Conformal Conic Projection
call GetData(Me%SP2, Me%ObjEnterData, flag2, &
keyword = 'SP2', &
default = 60., &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR130'
#ifdef _USE_PROJ4
#else
if (flag1 + flag2 == 2) then !estao la
Me%UseLambert = .TRUE.
else
Me%UseLambert = .FALSE.
endif
#endif
!Reads ZoneLong
if (Me%CoordType == UTM_ ) then
call ComputeGridZone (dble(Me%Longitude), dble(Me%Latitude), Me%Grid_Zone)
Me%ZoneLong = Me%Grid_Zone(1)
Me%ZoneLat = Me%Grid_Zone(2)
call GetData(ZoneLong, Me%ObjEnterData, flag, &
keyword = 'ZONE', &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR140'
if (flag == 1) then
if (ZoneLong /= Me%ZoneLong) stop 'ConstructGlobalVariables - HorizontalGrid - ERR150'
endif
else if (Me%CoordType .EQ. MIL_PORT_) then
Me%ZoneLong = PORTUGUESE_UTM_ZONE_
else
Me%ZoneLong = null_int
endif
!Reads ORIGIN
call GetData(AuxReal,Me%ObjEnterData, flag, &
keyword = 'ORIGIN', &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR160'
if (flag /= 2 ) stop 'ConstructGlobalVariables - HorizontalGrid - ERR170'
Me%Xorig = AuxReal(1)
Me%Yorig = AuxReal(2)
!Reads GridAngle
call GetData(Me%Grid_Angle, Me%ObjEnterData, flag, &
keyword = 'GRID_ANGLE', &
default = 0.0, &
ClientModule = 'HorizontalGrid', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR180'
!Allocates XX and YY
allocate(Me%XX(Me%Size%JLB+1 : Me%Size%JUB+1))
allocate(Me%YY(Me%Size%ILB+1 : Me%Size%IUB+1))
allocate(Me%XX_IE(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%YY_IE(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%LatitudeConn(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%LongitudeConn(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%DefineCellsMap(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%XX_AlongGrid(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%YY_AlongGrid(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
Me%XX = null_real
Me%YY = null_real
Me%XX_IE = null_real
Me%YY_IE = null_real
Me%LatitudeConn = null_real
Me%LongitudeConn = null_real
Me%DefineCellsMap(:,:) = 1
!Reads XX_YE, YY_IE
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR190'
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
BeginCornersXY, EndCornersXY, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR200'
BF: if (BlockFound) then
Me%CornersXYInput = .true.
if (Me%CoordType == CIRCULAR_ .and. Me%CoordType == GEOG_) then
stop 'ConstructGlobalVariables - HorizontalGrid - ERR210'
endif
line = FirstLine
do i = Me%GlobalWorkSize%ILB, Me%GlobalWorkSize%IUB+1
do j = Me%GlobalWorkSize%JLB, Me%GlobalWorkSize%JUB+1
line = line+1
!Last line found before end?
if (line == LastLine) then
write(*,*)
write(*,*) 'Error reading Corners X Y'
stop 'ConstructGlobalVariables - HorizontalGrid - ERR220'
end if
!Reads XX_IE , YY_IE
call GetData(AuxReal, Me%ObjEnterData, flag, &
Buffer_Line = Line, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR230'
if (Me%DDecomp%MasterOrSlave) then
if (i>= Me%DDecomp%HaloMap%ILB .and. &
i<= Me%DDecomp%HaloMap%IUB+1) then
ii = i - Me%DDecomp%HaloMap%ILB + 1
else
cycle
endif
else
ii = i
endif
if (Me%DDecomp%MasterOrSlave) then
if (j>= Me%DDecomp%HaloMap%JLB .and. &
j<= Me%DDecomp%HaloMap%JUB+1) then
jj = j - Me%DDecomp%HaloMap%JLB + 1
else
cycle
endif
else
jj = j
endif
Me%XX_IE(ii, jj) = AuxReal(1)
Me%YY_IE(ii, jj) = AuxReal(2)
end do
end do
if (Me%CoordType == CIRCULAR_ .or. Me%CornersXYInput) then
Me%Distortion = .true.
Me%RegularRotation = .false.
else
Me%Distortion = .false.
if (abs(Me%Grid_Angle)>0.) then
Me%RegularRotation = .true.
else
Me%RegularRotation = .false.
endif
endif
allocate(Me%DefineFacesUMap(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%DefineFacesVMap(Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
allocate(Me%DefineCrossMap (Me%Size%ILB : Me%Size%IUB, Me%Size%JLB : Me%Size%JUB))
Me%DefineFacesUMap(:,:) = 1
Me%DefineFacesVMap(:,:) = 1
Me%DefineCrossMap (:,:) = 1
Me%NotDefinedCells = .false.
!If one corner has FillValue Real, the point is not to consider
do i = Me%WorkSize%ILB, Me%WorkSize%IUB
do j = Me%WorkSize%JLB, Me%WorkSize%JUB
if (Me%XX_IE(i , j ) <= FillValueReal / 2. .or. &
Me%XX_IE(i , j+1) <= FillValueReal / 2. .or. &
Me%XX_IE(i+1, j+1) <= FillValueReal / 2. .or. &
Me%XX_IE(i+1, j ) <= FillValueReal / 2.) then
Me%DefineCellsMap(i, j) = 0
Me%NotDefinedCells = .true.
endif
end do
end do
do i = Me%WorkSize%ILB , Me%WorkSize%IUB
if (Me%DefineCellsMap(i, Me%WorkSize%JLB) == 0) Me%DefineFacesUMap(i, Me%WorkSize%JLB ) = 0
if (Me%DefineCellsMap(i, Me%WorkSize%JUB) == 0) Me%DefineFacesUMap(i, Me%WorkSize%JUB+1) = 0
do j = Me%WorkSize%JLB+1, Me%WorkSize%JUB
!if (Me%XX_IE(i , j ) <= FillValueReal / 2. .or. &
! Me%XX_IE(i+1, j ) <= FillValueReal / 2. ) then
! Me%DefineFacesUMap(i,j) = 0
!endif
if (Me%DefineCellsMap(i, j-1) == 0 .and. Me%DefineCellsMap(i, j) == 0) Me%DefineFacesUMap(i,j) = 0
end do
end do
do j = Me%WorkSize%JLB, Me%WorkSize%JUB
if (Me%DefineCellsMap(Me%WorkSize%ILB, j) == 0) Me%DefineFacesVMap(Me%WorkSize%ILB,j ) = 0
if (Me%DefineCellsMap(Me%WorkSize%IUB, j) == 0) Me%DefineFacesVMap(Me%WorkSize%IUB+1,j) = 0
do i = Me%WorkSize%ILB+1, Me%WorkSize%IUB
!if (Me%XX_IE(i , j+1) <= FillValueReal / 2. .or. &
! Me%XX_IE(i , j ) <= FillValueReal / 2. ) then
! Me%DefineFacesVMap(i,j) = 0
!endif
if (Me%DefineCellsMap(i-1, j) == 0 .and. Me%DefineCellsMap(i, j) == 0) Me%DefineFacesVMap(i,j) = 0
end do
end do
do i = Me%WorkSize%ILB+1, Me%WorkSize%IUB
do j = Me%WorkSize%JLB+1, Me%WorkSize%JUB
!if (Me%XX_IE(i , j ) <= FillValueReal / 2. ) then
! Me%DefineCrossMap(i, j) = 0
!endif
if (Me%DefineCellsMap(i-1, j ) == 0 .and. Me%DefineCellsMap(i, j ) == 0 .or. &
Me%DefineCellsMap(i-1, j-1) == 0 .and. Me%DefineCellsMap(i, j-1) == 0) then
Me%DefineCrossMap(i, j) = 0
endif
end do
end do
call Block_Unlock(Me%ObjEnterData, ClientNumber, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR240'
else BF
if (abs(Me%Grid_Angle)>0.) then
Me%RegularRotation = .true.
else
Me%RegularRotation = .false.
endif
!Check if the spacing in Y is constant
call GetData(ConstantSpacingY, &
Me%ObjEnterData, flag, &
SearchType = FromFile, &
keyword ='CONSTANT_SPACING_Y', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR250'
iconY: if (ConstantSpacingY) Then
!Get grid origin
call GetData(DY, &
Me%ObjEnterData, flag, &
SearchType = FromFile, &
keyword ='DY', &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR260'
XY_Aux = -DY
ii = 0
do i = Me%GlobalWorkSize%ILB, Me%GlobalWorkSize%IUB + 1
XY_Aux = XY_Aux + DY
if (Me%DDecomp%MasterOrSlave) then
if (i>= Me%DDecomp%HaloMap%ILB .and. &
i<= Me%DDecomp%HaloMap%IUB+1) then
ii = ii + 1
else
cycle
endif
else
ii = i
endif
Me%YY(ii) = XY_Aux
end do
else iconY
!Reads YY
call RewindBuffer(Me%ObjEnterData, STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR270'
call ExtractBlockFromBuffer(Me%ObjEnterData, ClientNumber, &
BeginYY, EndYY, BlockFound, &
FirstLine = FirstLine, LastLine = LastLine, &
STAT = STAT_CALL)
if (STAT_CALL /= SUCCESS_) stop 'ConstructGlobalVariables - HorizontalGrid - ERR280'
if (BlockFound) then
line = FirstLine
ii = 0
do i = Me%GlobalWorkSize%ILB, Me%GlobalWorkSize%IUB+1
line = line+1
!Last line found before end?
if (line == LastLine) then
write