Skip to content

Commit

Permalink
Merge pull request #12 from Mohid-Water-Modelling-System/master
Browse files Browse the repository at this point in the history 
update my local repository
  • Loading branch information
@sofiasaraiva
sofiasaraiva authored May 10, 2021
2 parents df507bb + 780efe6 commit 829b8a1
Show file tree
Hide file tree
Showing 26 changed files with 2,781 additions and 297 deletions.
308 changes: 293 additions & 15 deletions Software/ConvertToHDF5/ModuleNetCDFCF_2_HDF5MOHID.F90
View file

Large diffs are not rendered by default.

26 changes: 22 additions & 4 deletions Software/MOHIDBase1/ModuleDrawing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -969,17 +969,25 @@ subroutine SetLimitsPolygon(Polygon)

d1: do i=1, Polygon%Count

if (Polygon%VerticesF(i)%X < Polygon%Limits%Left) &
if (Polygon%VerticesF(i)%X < Polygon%Limits%Left .and. &
Polygon%VerticesF(i)%X > FillValueReal / 2.) then
Polygon%Limits%Left = Polygon%VerticesF(i)%X
endif

if (Polygon%VerticesF(i)%X > Polygon%Limits%Right) &
if (Polygon%VerticesF(i)%X > Polygon%Limits%Right .and. &
Polygon%VerticesF(i)%X > FillValueReal / 2.) then
Polygon%Limits%Right = Polygon%VerticesF(i)%X
endif

if (Polygon%VerticesF(i)%Y < Polygon%Limits%Bottom) &
if (Polygon%VerticesF(i)%Y < Polygon%Limits%Bottom .and. &
Polygon%VerticesF(i)%Y > FillValueReal / 2.) then
Polygon%Limits%Bottom = Polygon%VerticesF(i)%Y
endif

if (Polygon%VerticesF(i)%Y > Polygon%Limits%Top) &
if (Polygon%VerticesF(i)%Y > Polygon%Limits%Top .and. &
Polygon%VerticesF(i)%Y > FillValueReal / 2.) then
Polygon%Limits%Top = Polygon%VerticesF(i)%Y
endif

enddo d1

Expand Down Expand Up @@ -1504,6 +1512,16 @@ logical function IsPointInsidePolygon(Point, Polygon)

do1: do CurrentVertix = 1, Polygon%Count - 1

if (Polygon%VerticesF(CurrentVertix)%X == FillValueReal) then
write (*,*) 'Vertice with invalid X coordinate'
stop 'IsPointInsidePolygon - ModuleDrawing - ERR10'
endif

if (Polygon%VerticesF(CurrentVertix)%Y == FillValueReal) then
write (*,*) 'Vertice with invalid Y coordinate'
stop 'IsPointInsidePolygon - ModuleDrawing - ERR20'
endif

!construct segment
Segment%StartAt%X = Polygon%VerticesF(CurrentVertix)%X
Segment%StartAt%Y = Polygon%VerticesF(CurrentVertix)%Y
Expand Down
132 changes: 106 additions & 26 deletions Software/MOHIDBase1/ModuleFunctions.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,6 +249,7 @@ Module ModuleFunctions

!Polygon
public :: RelativePosition4VertPolygon
public :: RelativePosition4VertPolygonIWD
public :: PolygonArea
public :: FromGeo2Meters

Expand Down Expand Up @@ -7928,95 +7929,172 @@ subroutine RelativePosition4VertPolygon(Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye,

!Local-------------------------------------------------------------------
real(8) :: DXdc, DYac, DXba, DYbd, DXef, DYeg
real(8) :: MinDx, SumAux
real(8) :: MinDx
real(8) :: a1, b1, a2, b2, a3, b3, a4, b4
real(8) :: Seg_ac, Seg_dc,Seg_hc, Seg_ic
real(8) :: Xf, Yf, Xg, Yg, Xh, Yh, Xi, Yi, TgX, TgY
real(8) :: SumAux
real(8) :: XaR8, YaR8, XbR8, YbR8, XcR8, YcR8, XdR8, YdR8, XeR8, YeR8
!Begin-------------------------------------------------------------------

XaR8 = dble(Xa); YaR8 = dble(Ya); XbR8 = dble(Xb); YbR8 = dble(Yb); XcR8 = dble(Xc)
YcR8 = dble(Yc); XdR8 = dble(Xd); YdR8 = dble(Yd); XeR8 = dble(Xe); YeR8 = dble(Ye)




!the four segments of the cell
DXdc = XdR8 - XcR8
DXba = XbR8 - XaR8
DYac = YaR8 - YcR8
DYbd = YbR8 - YdR8



SumAux = abs(DXdc + DXba + DYac + DYbd)
MinDx = SumAux * 1.e-16


if (SumAux > 0) then
MinDx = SumAux * 1.e-16
else
MinDx = 1.e-16
endif

if (abs(DXdc)<MinDx) DXdc = MinDx
if (abs(DXba)<MinDx) DXba = MinDx
if (abs(DYac)<MinDx) DYac = MinDx
if (abs(DYbd)<MinDx) DYbd = MinDx

a1 = (XdR8*YcR8 - XcR8*YdR8) / DXdc
b1 = ( YdR8 - YcR8) / DXdc

a2 = (XbR8*YaR8 - XaR8*YbR8) / DXba
b2 = ( YbR8 - YaR8) / DXba

a3 = (XcR8*YaR8 - XaR8*YcR8) / DYac
b3 = ( XaR8 - XcR8) / DYac

a4 = (XdR8*YbR8 - XbR8*YdR8) / DYbd
b4 = ( XbR8 - XdR8) / DYbd

!intersection points
if (b2/=b1) then

!F point intersection point of X faces
Xf = (a1 - a2) / (b2 - b1)
Yf = a1 + b1*Xf

!H point intersection with segment CA
DXef = XeR8 - Xf

if (abs(DXef)==0.) DXef = MinDx
TgX = (YeR8 - Yf) / DXef
else
!H point intersection with segment CA
TgX = b1
endif

Yh = (YeR8 + (a3 - XeR8) * TgX) / (1. - b3*TgX)
Xh = a3 + b3 * Yh


if (b4/=b3) then

!G point intersection point of X faces
Yg = (a3 - a4) / (b4 - b3)
Xg = a3 + b3*Yg

!i point intersection with segment CA
DYeg = YeR8 - Yg

if (abs(DYeg)==0.) DYeg = MinDx
TgY = (XeR8 - Xg) / DYeg
else

!i point intersection with segment CA
TgY = b3
endif

Xi = (XeR8 + (a1 - YeR8) * TgY) / (1. - b1*TgY)
Yi = a1 + b1 * Xi


Seg_ac = sqrt((XaR8-XcR8)*(XaR8-XcR8) + (YaR8-YcR8)*(YaR8-YcR8))
Seg_dc = sqrt((XdR8-XcR8)*(XdR8-XcR8) + (YdR8-YcR8)*(YdR8-YcR8))
Seg_hc = sqrt((Xh -XcR8)*(Xh -XcR8) + (Yh -YcR8)*(Yh -YcR8))
Seg_ic = sqrt((Xi -XcR8)*(Xi -XcR8) + (Yi -YcR8)*(Yi -YcR8))

Xex = Seg_ic / Seg_dc
Yey = Seg_hc / Seg_ac

if (Xex < 0. .or. Xex > 1 .or. Yey < 0. .or. Yey > 1) then
call RelativePosition4VertPolygonIWD(Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye, Xex, Yey)
endif


end subroutine RelativePosition4VertPolygon
!--------------------------------------------------------------------------

subroutine RelativePosition4VertPolygonIWD(Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye, Xex, Yey)

!Arguments---------------------------------------------------------------
real :: Xa, Ya, Xb, Yb, Xc, Yc, Xd, Yd, Xe, Ye, Xex, Yey

!Local-------------------------------------------------------------------
real(8) :: SumAux
real(8) :: XaR8, YaR8, XbR8, YbR8, XcR8, YcR8, XdR8, YdR8, XeR8, YeR8
real(8) :: da, db, dc, dd, pw
!Begin-------------------------------------------------------------------

XaR8 = dble(Xa); YaR8 = dble(Ya); XbR8 = dble(Xb); YbR8 = dble(Yb); XcR8 = dble(Xc)
YcR8 = dble(Yc); XdR8 = dble(Xd); YdR8 = dble(Yd); XeR8 = dble(Xe); YeR8 = dble(Ye)

pw = 2

!using a IDW - Inverse distance weighting method
!(Xa,Ya) = (Xex = 0, Yey =1)
!(Xb,Yb) = (Xex = 1, Yey =1)
!(Xc,Yc) = (Xex = 0, Yey =0)
!(Xd,Yd) = (Xex = 1, Yey =0)
!
! a_______b
! | |
! | |
! c_______d

da = sqrt((XeR8-XaR8)**2+(YeR8-YaR8)**2)
db = sqrt((XeR8-XbR8)**2+(YeR8-YbR8)**2)
dc = sqrt((XeR8-XcR8)**2+(YeR8-YcR8)**2)
dd = sqrt((XeR8-XdR8)**2+(YeR8-YdR8)**2)

if (da == 0) then
Xex = 0
Yey = 1
elseif (db == 0) then
Xex = 1
Yey = 1
elseif (dc == 0) then
Xex = 0
Yey = 0
elseif (dd == 0) then
Xex = 1
Yey = 0
else
SumAux = ((1./da)**pw+(1./db)**pw+(1./dc)**pw+(1./dd)**pw)
Xex = ((1./db)**pw+(1./dd)**pw) / SumAux
Yey = ((1./da)**pw+(1./db)**pw) / SumAux
endif


if (Xex < 0. .or. Xex > 1) then
stop 'Function - RelativePosition4VertPolygonIWD - ERR10'
endif

if (Yey < 0. .or. Yey > 1) then
stop 'Function - RelativePosition4VertPolygonIWD - ERR20'
endif


end subroutine RelativePosition4VertPolygonIWD
!--------------------------------------------------------------------------



!! Public-domain function by Darel Rex Finley, 2006.
Expand Down Expand Up @@ -13746,7 +13824,8 @@ end subroutine DischargeFluxV
!>@param[in] Flow, DischargeConnection, VelFather, VelSon, AreaU, DecayTime, VelDT, CoefCold
subroutine Offline_DischargeFluxU(Flow, DischargeConnection, VelFather, VelSon, AreaU, DecayTime, VelDT, CoefCold)
!Arguments--------------------------------------------------------------------------
real, dimension(:, :, :), pointer, intent(IN) :: VelFather, VelSon, AreaU
real, dimension(:, :, :), pointer, intent(IN) :: VelFather, AreaU
real, dimension(:, :, :), allocatable, intent(IN) :: VelSon
real(8), dimension(:) , intent(INOUT) :: Flow
integer, dimension(:, :) , allocatable, intent(IN) :: DischargeConnection
real , dimension(:, :) , pointer, intent(IN) :: DecayTime
Expand Down Expand Up @@ -13801,7 +13880,8 @@ end subroutine Offline_DischargeFluxU
!>@param[in] Flow, DischargeConnection, VelFather, VelSon, AreaU, DecayTime, VelDT, CoefCold
subroutine Offline_DischargeFluxV(Flow, DischargeConnection, VelFather, VelSon, AreaV, DecayTime, VelDT, CoefCold)
!Arguments--------------------------------------------------------------------------
real, dimension(:, :, :), pointer, intent(IN) :: VelFather, VelSon, AreaV
real, dimension(:, :, :), pointer, intent(IN) :: VelFather, AreaV
real, dimension(:, :, :), allocatable, intent(IN) :: VelSon
real(8), dimension(:) , intent(OUT) :: Flow
integer, dimension(:, :) , allocatable, intent(IN) :: DischargeConnection
real , dimension(:, :) , pointer, intent(IN) :: DecayTime
Expand Down
11 changes: 8 additions & 3 deletions Software/MOHIDBase1/ModuleGlobalData.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,10 +76,12 @@ Module ModuleGlobalData
end interface SetError

!Parameter-----------------------------------------------------------------
integer, parameter :: MaxModules = 98
integer, parameter :: MaxModules = 100

#ifdef _INCREASE_MAXINSTANCES
#if defined(_INCREASE_MAXINSTANCES)
integer, parameter :: MaxInstances = 2000
#elif defined(_INCREASE_MAXINSTANCES_EXTRA)
integer, parameter :: MaxInstances = 20000
#else
integer, parameter :: MaxInstances = 500
#endif
Expand Down Expand Up @@ -1980,6 +1982,8 @@ Module ModuleGlobalData
integer, parameter :: mLitter_ = 96
integer, parameter :: mTwoWay_ = 97
integer, parameter :: mOutputGrid_ = 98
integer, parameter :: mMeshGlue_ = 99
integer, parameter :: mDelftFM_2_MOHID_ = 100

!Domain decomposition
integer, parameter :: WestSouth = 1
Expand Down Expand Up @@ -2108,7 +2112,8 @@ Module ModuleGlobalData
T_Module(mSediment_ , "Sediment" ), T_Module(mReservoirs_ , "Reservoirs" ), &
T_Module(mIrrigation_ , "Irrigation" ), T_Module(mTURBINE_ , "Turbine" ), &
T_Module(mLitter_ , "Litter" ), T_Module(mTwoWay_ , "TwoWay" ), &
T_Module(mOutputGrid_ , "OuputGrid" )/)
T_Module(mOutputGrid_ , "OutputGrid" ), T_Module(mMeshGlue_ , "MeshGlue" ), &
T_Module(mDelftFM_2_MOHID_ , "DelftFM_2_MOHID" )/)


!Variables
Expand Down
11 changes: 8 additions & 3 deletions Software/MOHIDBase1/ModuleHDF5.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7458,7 +7458,7 @@ recursive subroutine CheckAllDataSets (IDIn, GroupNameIn)

allocate (offset_in (rank))
allocate (count_in (rank))
offset_in(:) = 1
offset_in(:) = 0

count_in (:) = 1

Expand Down Expand Up @@ -7506,7 +7506,7 @@ recursive subroutine CheckAllDataSets (IDIn, GroupNameIn)

NumType = H5T_NATIVE_REAL

call h5dread_f (dset_id, NumType, ValueOut(1),&
call h5dread_f (dset_id, NumType, ValueOut(1:1),&
dims_mem, STAT, memspace_id, space_id)
if (STAT /= SUCCESS_) then
write(*,*) "FileName =",trim(Me%Filename)
Expand All @@ -7529,6 +7529,11 @@ recursive subroutine CheckAllDataSets (IDIn, GroupNameIn)
if (STAT /= SUCCESS_) then
stop 'CheckAllDataSets - ModuleHDF5 - ERR110'
endif

call h5sclose_f (memspace_id, STAT)
if (STAT /= SUCCESS_) then
stop 'CheckAllDataSets - ModuleHDF5 - ERR115'
endif



Expand Down Expand Up @@ -8119,7 +8124,7 @@ subroutine LocateObjHDF5 (HDF5ID)
enddo

if (.not. associated(Me)) then
write(*,*) Me%FileName
write(*,*) "HDF5ID =", HDF5ID
stop 'ModuleHDF5 - LocateObjHDF5 - ERR01'
endif

Expand Down
Loading

0 comments on commit 829b8a1

Please sign in to comment.