!----- AGPL --------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2015.
!
! This file is part of Delft3D (D-Flow Flexible Mesh component).
!
! Delft3D is free software: you can redistribute it and/or modify
! it under the terms of the GNU Affero General Public License as
! published by the Free Software Foundation version 3.
!
! Delft3D 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 Affero General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License
! along with Delft3D. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D",
! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting
! Deltares, and remain the property of Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: rest.F90 43422 2015-12-04 17:03:07Z pijl $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/rest.F90 $
!> Checks whether lines 1-2 and 3-4 intersect.
!! @param[in] x1,y1,x2,y2,x3,y3,x4,y4 x- and y-coords of line endpoints.
!! @param[out] jacros 1 if lines cross (intersect), 0 if not.
!! @param[out] sl lambda in [0,1] on line segment 1-2 (outside [0,1] if no intersection). Unchanged if no intersect!!
!! @param[out] sm lambda in [0,1] on line segment 3-4 (outside [0,1] if no intersection). Unchanged if no intersect!!
!! @param[out] xcr,ycr x-coord. of intersection point.
SUBROUTINE CROSS(x1, y1, x2, y2, x3, y3, x4, y4, JACROS,SL,SM,XCR,YCR,CRP)
use m_missing
implicit none
double precision, intent(inout) :: crp !< crp (in)==-1234 will make crp (out) non-dimensional
double precision :: det
double precision :: eps
integer :: jacros, jamakenondimensional
double precision :: sl
double precision :: sm
double precision, intent(in) :: x1, y1, x2, y2, x3, y3, x4, y4
double precision :: x21, y21, x31, y31, x43, y43, xcr, ycr
double precision, external :: getdx, getdy
! safety check on crp (in)
if ( isnan(crp) ) then
crp = 0d0
end if
! Set defaults for no crossing at all:
jamakenondimensional = 0
if ( abs(crp+1234d0).lt.0.5d0 ) then
jamakenondimensional = 1
crp = 0d0
endif
JACROS = 0
EPS = 0.00001d0
SL = DMISS
SM = DMISS
! SL = LABDA TUSSEN 0 EN 1 OP EERSTE PAAR
! Sm = LABDA TUSSEN 0 EN 1 OP TWEEDE PAAR
call getdxdy(x1,y1,x2,y2,x21,y21)
call getdxdy(x3,y3,x4,y4,x43,y43)
call getdxdy(x1,y1,x3,y3,x31,y31)
!X21 = getdx(x1,y1,x2,y2)
!Y21 = getdy(x1,y1,x2,y2)
!X43 = getdx(x3,y3,x4,y4)
!Y43 = getdy(x3,y3,x4,y4)
!X31 = getdx(x1,y1,x3,y3)
!Y31 = getdy(x1,y1,x3,y3)
DET = X43*Y21 - Y43*X21
! SPvdP: make eps have proper dimension
EPS = max(EPS*MAXVAL((/ X21,Y21,X43,Y43,X31,Y31 /)), tiny(0d0))
IF (ABS(DET) .LT. EPS) THEN
RETURN
ELSE
SM = (Y31*X21 - X31*Y21) / DET
IF (ABS(X21) .GT. EPS) THEN
SL = (SM*X43 + X31) / X21
ELSE IF (ABS(Y21) .GT. EPS) THEN
SL = (SM*Y43 + Y31) / Y21
ELSE
SL = 0d0
ENDIF
IF (SM .GE. 0d0 .AND. SM .LE. 1d0 .AND. &
SL .GE. 0d0 .AND. SL .LE. 1d0) THEN
JACROS = 1
ENDIF
XCR = X1 + SL*(X2-X1)
YCR = Y1 + SL*(Y2-Y1)
if ( jamakenondimensional.eq.1 ) then ! make crp non-dimensional (for spline2curvi)
CRP = -DET / ( sqrt(x21**2+y21**2) * sqrt(x43**2 + y43**2) + 1d-8 )
else
CRP = -DET
end if
ENDIF
RETURN
END
SUBROUTINE CROSSinbox (x1, y1, x2, y2, x3, y3, x4, y4, JACROS,SL,SM,XCR,YCR,CRP) ! only if overlap
use m_missing
implicit none
double precision, intent(inout) :: crp ! crp (in)==-1234 will make crp (out) non-dimensional
integer :: jacros
double precision, intent(in) :: x1, y1, x2, y2, x3, y3, x4, y4
double precision, intent(out) :: SL,SM,XCR,YCR
double precision :: x1min, x1max, y1min, y1max, x3min, x3max, y3min, y3max
! Set defaults for no crossing at all:
JACROS = 0
x1min = min(x1,x2); x1max = max(x1,x2)
y1min = min(y1,y2); y1max = max(y1,y2)
x3min = min(x3,x4); x3max = max(x3,x4)
y3min = min(y3,y4); y3max = max(y3,y4)
if (x1max < x3min) return
if (x1min > x3max) return
if (y1max < y3min) return
if (y1min > y3max) return
call CROSS (x1, y1, x2, y2, x3, y3, x4, y4, JACROS,SL,SM,XCR,YCR,CRP)
RETURN
END
SUBROUTINE dCROSS(X1,Y1,X2,Y2,X3,Y3,X4,Y4,JACROS,SL,SM,XCR,YCR,CRP) ! liggen 3 en 4 aan weerszijden van lijn 12
implicit none
double precision :: det
double precision :: eps
DOUBLE PRECISION :: X1,Y1,X2,Y2,X3,Y3,X4,Y4,SL,SM,XCR,YCR,CRP
INTEGER :: JACROS
DOUBLE PRECISION :: X21, Y21, X43, Y43, X31, Y31, getdx, getdy
JACROS = 0
EPS = 0.00001d0
! SL = LABDA TUSSEN 0 EN 1 OP EERSTE PAAR
! Sm = LABDA TUSSEN 0 EN 1 OP TWEEDE PAAR
call getdxdy(x1,y1,x2,y2,x21,y21)
call getdxdy(x3,y3,x4,y4,x43,y43)
call getdxdy(x1,y1,x3,y3,x31,y31)
!X21 = getdx(x1,y1,x2,y2)
!Y21 = getdy(x1,y1,x2,y2)
!X43 = getdx(x3,y3,x4,y4)
!Y43 = getdy(x3,y3,x4,y4)
!X31 = getdx(x1,y1,x3,y3)
!Y31 = getdy(x1,y1,x3,y3)
DET = X43*Y21 - Y43*X21
IF (ABS(DET) .LT. EPS) THEN
RETURN
ELSE
SM = (Y31*X21 - X31*Y21) / DET
IF (ABS(X21) .GT. EPS) THEN
SL = (SM*X43 + X31) / X21
ELSE IF (ABS(Y21) .GT. EPS) THEN
SL = (SM*Y43 + Y31) / Y21
ELSE
SL = 0d0
ENDIF
XCR = X1 + SL*(X2-X1)
YCR = Y1 + SL*(Y2-Y1)
CRP = -DET
IF (SM >= 0d0 .AND. SM <= 1d0) then
JACROS = 1
ENDIF
ENDIF
RETURN
END subroutine dcross
double precision FUNCTION DISLIN(X,Y,N,XX,YY,TV)
! AFSTAND VAN PUNT XX,YY TOT LIJN MET PARM TV
implicit none
integer :: n
double precision :: tv
double precision :: xv
double precision :: xx
double precision :: yv
double precision :: yy
double precision :: dbdistance
double precision :: X(N), Y(N)
TV = MAX(0d0,MIN(TV,N-1d0))
CALL LINT(X,Y,N,TV,XV,YV)
dislin = dbDISTANCE(XV,YV,XX,YY )
RETURN
END function dislin
SUBROUTINE LINT(X,Y,N,TV,XV,YV)
implicit none
integer :: n
integer :: n1
integer :: n2
integer :: ntv
double precision :: t
double precision :: tv
double precision :: xv
double precision :: yv
! Lineaire interpolatie op TV in lijn
double precision :: X(N), Y(N)
NTV = INT(TV)
T = TV - NTV
N1 = NTV + 1
N2 = N1 + 1
XV = (1-T)*X(N1) + T*X(N2)
YV = (1-T)*Y(N1) + T*Y(N2)
RETURN
END subroutine lint
SUBROUTINE GOLDLN(AX,BX,CX,TOL,XMIN,P,Q,N,XX,YY,DIS)
implicit none
double precision :: ax
double precision :: bx
double precision :: c
double precision :: cx
double precision :: dis
double precision :: f0
double precision :: f1
double precision :: f2
double precision :: f3
integer :: n
double precision :: r
double precision :: tol
double precision :: x0
double precision :: x1
double precision :: x2
double precision :: x3
double precision :: xmin
double precision :: xx
double precision :: yy
PARAMETER (R=.61803399,C=.38196602)
! EENDIMENSIONAAL ZOEKEN VAN 'GEBRACKED' MINIMUM
double precision :: P(N), Q(N)
double precision :: dislin
X0=AX
X3=CX
IF(ABS(CX-BX).GT.ABS(BX-AX))THEN
X1=BX
X2=BX+C*(CX-BX)
ELSE
X2=BX
X1=BX-C*(BX-AX)
ENDIF
! F1=F(X1)
! F1=DIST(P,P2,Q,Q2,XX,YY,X1,N)
F1=DISLIN(P,Q,N,XX,YY,X1)
! F2=F(X2)
! F2=DIST(P,P2,Q,Q2,XX,YY,X2,N)
F2=DISLIN(P,Q,N,XX,YY,X2)
1 IF(ABS(X3-X0).GT.TOL*(ABS(X1)+ABS(X2)))THEN
IF(F2.LT.F1)THEN
X0=X1
X1=X2
X2=R*X1+C*X3
F0=F1
F1=F2
! F2=DIST(P,P2,Q,Q2,XX,YY,X2,N)
F2=DISLIN(P,Q,N,XX,YY,X2)
! F2=F(X2)
ELSE
X3=X2
X2=X1
X1=R*X2+C*X0
F3=F2
F2=F1
! F1=F(X1)
F1=DISLIN(P,Q,N,XX,YY,X1)
! F1=DIST(P,P2,Q,Q2,XX,YY,X1,N)
ENDIF
GOTO 1
ENDIF
IF(F1.LT.F2)THEN
DIS =F1
XMIN=X1
ELSE
DIS =F2
XMIN=X2
ENDIF
RETURN
END
SUBROUTINE INCREASELAN(N)
USE M_MISSING
use unstruc_messages
use m_alloc
use m_landboundary
integer :: n
integer :: ierr
IF (N < MAXLAN) RETURN
MAXLAN = MAX(50000,INT(1.2d0*N))
call realloc(xlan, MAXLAN, stat=ierr, fill=dxymis)
CALL AERR('xlan(maxlan)', IERR, maxlan)
call realloc(ylan, MAXLAN, stat=ierr, fill=dxymis)
CALL AERR('ylan(maxlan)', IERR, maxlan)
call realloc(zlan, MAXLAN, stat=ierr, fill=dxymis)
CALL AERR('zlan(maxlan)', IERR, maxlan)
call realloc(nclan, MAXLAN, stat=ierr, fill=0)
CALL AERR('nclan(maxlan)', IERR, maxlan/2)
END subroutine increaselan
SUBROUTINE REALAN( MLAN)
use m_polygon
use M_landboundary
USE M_MISSING
implicit none
integer :: i
integer :: mlan
integer :: ncl
integer :: newlin
integer :: nkol
integer :: nrow
integer :: ntot, n, k, kd, ku
double precision :: xlr
CHARACTER CHARMC*5, MATR*4, REC*132
DOUBLE PRECISION :: XL, YL, ZL
NTOT = 0
call increaselan(10000)
CALL READYY('READING land boundary',0d0)
10 CONTINUE
READ(MLAN,'(A)',END=777,ERR=887) MATR
IF (MATR(1:1) .EQ. '*') GOTO 10
READ(MLAN,'(A)',END = 777) REC
READ(REC,*,ERR = 666) NROW, NKOL
NEWLIN = 0
DO 20 I = 1,NROW
IF (NTOT .GE. MAXLAN-1) THEN
call increaselan(NTOT+1)
ENDIF
READ(MLAN,'(A)',END = 999) REC
NCL = 0
ZL = 0
if (nkol == 2) then
READ (REC,*,ERR=881) XL,YL
else if (nkol == 3) then
READ (REC,*,ERR=881) XL,YL,NCL
else if (nkol == 4) then
READ (REC,*,ERR=881) XL,YL,ZL,NCL
endif
XLR = XL
881 IF (XL .EQ. 999.999d0 .OR. XLR == 999.999d0) THEN
XL = dmiss
YL = dmiss
ZL = dmiss
NCL = 0
ENDIF
IF (NTOT == 0) THEN
NTOT = NTOT + 1
MXLAN = NTOT
XLAN(NTOT) = XL
YLAN(NTOT) = YL
ZLAN(NTOT) = ZL
NCLAN(NTOT) = NCL
ELSE IF (XL .ne. XLAN(NTOT) .or. YL .ne. YLAN(NTOT) ) THEN
NTOT = NTOT + 1
MXLAN = NTOT
XLAN(NTOT) = XL
YLAN(NTOT) = YL
ZLAN(NTOT) = ZL
NCLAN(NTOT) = NCL
ENDIF
IF (MOD(I,1000) .EQ. 0) THEN
CALL READYY(' ',MIN( 1d0,dble(I)/MAXLAN ) )
ENDIF
20 CONTINUE
NTOT = NTOT + 1
MXLAN = NTOT
XLAN(NTOT) = dmiss
YLAN(NTOT) = dmiss
ZLAN(NTOT) = dmiss
GOTO 10
777 CONTINUE
MXLAN = NTOT
CALL READYY(' ', 1d0)
CALL READYY(' ',-1d0)
call doclose (MLAN)
return
n = 1 ! remove double points in lineseg oriented files
xpl(n) = xlan(1) ; ypl(n) = ylan(1)
do k = 2,mxlan-1
kd = k - 1; ku = k + 1
if (xlan(k) == dmiss .and. xlan(kd) == xlan(ku) .and. ylan(kd) == ylan(ku) ) then
else
n = n + 1
xpl(n) = xlan(k) ; ypl(n) = ylan(k)
endif
enddo
n = n + 1
xpl(n) = xlan(mxlan) ; ypl(n) = ylan(mxlan)
npl = n
RETURN
666 CALL QNREADERROR('SEARCHING NROWS,NCOLS, BUT GETTING', REC, MLAN)
MXLAN = NTOT
CALL READYY(' ', 1d0)
CALL READYY(' ',-1d0)
call doclose (MLAN)
RETURN
888 CALL QNREADERROR('SEARCHING COORDINATES, BUT GETTING', REC, MLAN)
MXLAN = NTOT
CALL READYY(' ', 1d0)
CALL READYY(' ',-1d0)
call doclose (MLAN)
RETURN
887 CALL QNREADERROR('EXPECTING 4 CHAR, BUT GETTING', MATR, MLAN)
MXLAN = NTOT
CALL READYY(' ', 1d0)
CALL READYY(' ',-1d0)
call doclose (MLAN)
RETURN
999 CALL QNEOFERROR(MLAN)
MXLAN = NTOT
CALL READYY(' ', 1d0)
CALL READYY(' ',-1d0)
call doclose (MLAN)
RETURN
END
!> Read land boundary from world vector shoreline files (NetCDF format)
!! Global Self-consistent Hierarchical High-resolution Shorelines (GSHHS)
!! http://opendap.deltares.nl/thredds/catalog/opendap/noaa/gshhs/catalog.html
!! Directly stored in m_landboundary module variables.
subroutine read_land_boundary_netcdf(filename)
use M_landboundary
USE M_MISSING
use netcdf
implicit none
character(len=*), intent(in) :: fileName
double precision, dimension(:), allocatable :: x_lan, y_lan
integer, dimension(:), allocatable :: k
integer :: i
integer :: status, istart, istop
integer :: id_nc, id_lon, id_lat,id_npts, id_sep, id_k
integer :: npts, nsep
integer :: max_vertex, n_vertex
logical :: succes
! write(msgtxt,'("Reading file: ",a)') trim(get_basename( get_filename(filename) ) )
! Read the land boundary geometry from netcdf-file.
succes = .false.
status = nf90_open(trim(filename), NF90_NOWRITE, ncid=id_nc)
! if (status .ne. nf90_noerr) msgtxt = nf90_strerror(status)
status = nf90_inq_dimid(id_nc, "npoints", id_npts)
status = nf90_inquire_dimension(id_nc, id_npts, len=npts)
status = nf90_inq_dimid(id_nc, "segment_separators", id_sep)
status = nf90_inquire_dimension(id_nc, id_sep, len=nsep)
!call setTotalSteps(pbar, 5+nsep)
allocate(k(nsep))
status = nf90_inq_varid(id_nc, "k", id_k)
status = nf90_get_var(id_nc, id_k, k, count=(/ nsep /))
status = nf90_inq_varid(id_nc, "lon", id_lon)
! if (status .ne. nf90_noerr) msgtxt = nf90_strerror(status)
status = nf90_inq_varid(id_nc, "lat", id_lat)
! if (status .ne. nf90_noerr) msgtxt = nf90_strerror(status)
! Determine largest landboundary segment
max_vertex = 0
do i = 1, nsep-1
max_vertex = max(max_vertex, k(i+1)-k(i)-1 )
enddo
!allocate(x_lan(npts))
!allocate(y_lan(npts))
call increaselan(npts)
!call setProgress(pbar, 5)
status = nf90_get_var(id_nc, id_lon, xlan, count=(/ npts /))
status = nf90_get_var(id_nc, id_lat, ylan, count=(/ npts /))
!call multiFeatureSetCapacity(polylines, nsep-1+10)
do i = 1, nsep
! if (mod(i,100)==0) call setProgress(pbar, 5+i)
!istart = k(i ) + 1
!istop = k(i+1) - 1
!n_vertex = istop -istart + 1
! Replace NetCDF NaNs by our dmiss vals at segment separator positions, that's all.
xlan(k(i)) = dmiss
ylan(k(i)) = dmiss
!
!polyline = newpolyline()
!call addPointArray( polyline, n_vertex, x_lan(istart:istop), y_lan(istart:istop))
!
!call multifeatureaddfeature(polylines, polyline)
enddo
MXLAN = npts
status = nf90_close(id_nc)
!call setProgress(pbar, nsep)
!call free(pbar)
deallocate(k)
!deallocate(x_lan)
!deallocate(y_lan)
if (status==0) then
succes = .true.
endif
end subroutine read_land_boundary_netcdf
!> Read polygon file (or cross section/pli file) and store in global polygon.
!! File should contain Tekal block(s) with two or three columns.
!! The block names may be used for cross sections.
!! A dmiss line starts a new polyline without a name. Multiple dmiss lines are skipped.
SUBROUTINE REAPOL(MPOL, jadoorladen)
USE M_POLYGON
use network_data, only: netstat, NETSTAT_CELLS_DIRTY
USE M_MISSING
use m_alloc
use unstruc_messages
use unstruc_files
implicit none
integer :: mpol
integer, intent(in) :: jadoorladen !< Append to existing polygons (intended to read multiple crs files)
integer :: i, ipli, janampl
integer :: nkol
integer :: nrow
integer :: nmiss
double precision :: xx, yy, zz, dz1, dz2
CHARACTER(len=5) :: CHARMC
character(len=64) :: MATR
character(len=132) :: REC
janampl = 0
if (index (filenames(mpol),'crs' ) > 0) then
janampl = 1
endif
if (index (filenames(mpol),'CRS' ) > 0) then
janampl = 1
endif
ipli = 0
if (jadoorladen /= 1) then
if (.not. allocated(XPL)) allocate(XPL(1), YPL(1), ZPL(1))
XPL = XYMIS
YPL = XYMIS
ZPL = XYMIS
NPL = 0
if (allocated(nampli)) nampli = ' '
else
! Count number of existing polylines, to start storing new names at proper index.
do i=1,npl-1
if (xpl(i) == dmiss) then
ipli = ipli + 1
end if
end do
end if
call realloc(nampli,20, fill = ' ')
CALL READYY('READING POLYGON / land boundary / CRS-FILE',0d0)
10 CONTINUE
READ(MPOL,'(A)',END=999,ERR=888) MATR
IF (MATR(1:1) .EQ. '*' .or. len_trim(matr) == 0) GOTO 10
READ(MPOL,'(A)',END = 999) REC
READ(REC,*,ERR=888) NROW, NKOL
jaKol45 = 0
if (nkol < 2) then
CALL QNERROR('File should contain at least 2 or 3 columns, but got:', ' ', ' ') ! nkol)
goto 999
else if (nkol > 3) then
jaKol45 = 1
end if
CALL INCREASEPOL(NPL + NROW + 1, 1) ! previous pols (if any) + 1 dmiss + new polyline
11 ipli = ipli + 1 ! Start reading a new polyline
if (janampl == 1) then ! hk: only store names when called from reacrosssections
if (len_trim(matr) > 0) then
call realloc(nampli, int(1.2*ipli) + 1, fill = ' ')
nampli(ipli) = matr ! Temporarily store cross section name with polyline
endif
endif
if (npl > 0 .and. nrow > 0) then
! Separator element for subsequent named polylines
npl = npl + 1
xpl(npl) = dmiss
ypl(npl) = dmiss
zpl(npl) = dmiss
end if
I = 0
row: DO
I = I+1
if (I > NROW) exit
nmiss = 0
do ! Read a single line, or multiple until a NON-dmiss line is found
READ(MPOL,'(A)',END = 999) REC
ZZ = DMISS ; dz1 = dmiss; dz2 = dmiss
if (nkol >= 5) then
READ(REC,*,ERR=777) XX,YY,ZZ,dz1,dz2
else if (nkol == 4) then
READ(REC,*,ERR=777) XX,YY,ZZ,dz1
else if (nkol == 3) then
READ(REC,*,ERR=777) XX,YY,ZZ
else
READ(REC,*,ERR=777) XX,YY
end if
IF (XX .NE. dmiss .AND. XX .NE. 999.999d0) exit
nmiss = nmiss + 1
I = I+1
if (I > NROW) exit row ! Last row was also dmiss, now exit outer 'row' loop.
end do
if (nmiss > 0) then
backspace(MPOL) ! Last one was NON-dmiss, preceded by one or more dmiss lines
NROW = NROW-I+1 ! so reread this in the next polyline loop
MATR = ' ' ! Fake new Tekal block by decrementing NROW and dummy name in MATR
goto 11
else
NPL = NPL + 1
XPL(NPL) = XX
YPL(NPL) = YY
ZPL(NPL) = ZZ
if (jakol45 == 1) then
IF (.NOT. ALLOCATED(DZL) ) THEN
ALLOCATE ( DZL(MAXPOL), DZR(MAXPOL) )
ENDIF
DZL(NPL) = DZ1
DZR(NPL) = DZ2
endif
end if
IF (MOD(NPL,100) .EQ. 0) THEN
CALL READYY(' ',MIN( 1d0,dble(I)/MAXPOL ) )
ENDIF
end do row
GOTO 10
999 CONTINUE
! If last polyline in file had only dmisses, remove last separator element in xpl.
if (npl > 0) then
if (xpl(npl) == dmiss) then
npl = npl-1
end if
end if
CALL READYY(' ',-1d0)
call doclose (MPOL)
RETURN
888 CALL QNREADERROR('SEARCHING NROWS,NCOLS, BUT GETTING', REC, MPOL)
CALL READYY(' ',-1d0)
call doclose (MPOL)
RETURN
777 CALL QNREADERROR('READING COORDINATES BUT GETTING',REC,MPOL)
CALL READYY(' ',-1d0)
call doclose (MPOL)
RETURN
END subroutine reapol
SUBROUTINE WRIPOL(MPOL)
USE M_POLYGON
use m_missing
implicit none
integer :: mpol, numnampli
integer :: NCLAN(0)
double precision :: ZSH(0)
if (NPL<=0) return
numnampli = size(nampli)
if (zpl(1) == dmiss) then ! No third column for z-values
CALL WRILDB(MPOL, XPL, YPL, NPL, NCLAN, 0, ZSH, 0, nampli, 64, numnampli)
else
CALL WRILDB(MPOL, XPL, YPL, NPL, NCLAN, 0, ZPL, NPL, nampli, 64, numnampli)
end if
END subroutine wripol
SUBROUTINE WRILAN(MPOL)
USE M_LANDBOUNDARY
implicit none
integer :: mpol
integer :: mx
double precision, ALLOCATABLE :: XL(:), YL(:)
double precision :: ZL(0) ! no z-values
character(len=1) :: names(1) ! no names
MX = MAXLAN
ALLOCATE ( XL(MX), YL(MX))
XL (1:MXLAN) = XLAN(1:MXLAN)
YL (1:MXLAN) = YLAN(1:MXLAN)
names = ' '
CALL WRILDB(MPOL, XL, YL, MXLAN, nclan, MXLAN, ZL, 0, names, 1, 1)
DEALLOCATE (XL, YL)
END
!> Writes active cross sections to a polyline file.
subroutine wricrs(mpol)
use m_crosssections
use m_polygon
use m_missing
implicit none
integer :: mpol,i
call savepol()
call copycrosssectionstopol()
! npl = 0 ! Write traced polygons instead of original plis
! do i=1,ncrs
! xpl(npl+1:npl+crs(i)%len+1)=crs(i)%xk(1:crs(i)%len+1)
! ypl(npl+1:npl+crs(i)%len+1)=crs(i)%yk(1:crs(i)%len+1)
! npl = npl+crs(i)%len+2
! xpl(npl) = dmiss
! ypl(npl) = dmiss
! end do
! if (ncrs>0) npl = npl - 1 ! remove last separator
call wripol(mpol)
call restorepol()
end subroutine wricrs
!> Writes a polygon/land boundary/cross section file.
!! The polyline(s) are written as a sequence of Tekal blocks.
!! The name for each Tekal block can be specified, or is auto-generated
!! as 'L00x' otherwise.
SUBROUTINE WRILDB(MPOL, XSH, YSH, NSH, NCLAN, nnclan, ZSH, nzsh, names, namlen, nnam)
USE M_MISSING
implicit none
integer, intent(inout) :: mpol !< Open file pointer where to write to.
double precision, intent(in) :: XSH(NSH), YSH(NSH) !< Coordinates, polylines can be separated by dmiss value.
integer, intent(in) :: nsh !< Number of points in polyline.
character(len=namlen), intent(in) :: names(nnam) !< Names of all polylines, header of each Tekal Block.
integer, intent(in) :: namlen !< string length of names.
integer, intent(in) :: nnam !< Number of polyline names.
integer, intent(in) :: NCLAN(*) !< Third integer value for each point in XSH, optional: use nnclan=0 to ignore
integer, intent(in) :: nnclan !< Size of NCLAN, use 0 to ignore.
double precision, intent(in) :: ZSH(*) !< Third double value for each point in XSH, optional: use nzsh=0 to ignore
integer, intent(in) :: nzsh !< Size of ZSH, use 0 to ignore.
integer :: L
integer :: i, ipli, npli, mbna, ncol
integer, allocatable :: istart(:), iend(:)
character(len=max(namlen,10)) :: name
character(len=1) :: cdigits
character(len=40) :: rec
logical :: jaNCLAN, jaZSH
! Only include third column when size is equal to XSH array (or larger).
jaNCLAN = nNCLAN >= NSH
jaZSH = nZSH >= NSH
CALL READYY('Writing Polygon / Land Boundary FILE',0d0)
MBNA = 0
IF (MBNA > 0) call newfil(mbna, 'bna.bna')
if (NSH <= 0) goto 11
allocate(istart(nsh), iend(nsh))
! First, find starts and ends of all polylines (separated by dmiss line(s))
! such that each can be written as a named Tekal block later.
ipli = 0
i = 0
pli: do
i = i+1
if (i > nsh) exit pli
if (xsh(i) == dmiss) cycle pli
! Start of a new polyline found
ipli= ipli + 1
istart(ipli) = i
pts: do
i = i+1
if (i > nsh) exit pts
if (xsh(i) == dmiss) exit pts
end do pts
iend(ipli) = i-1
end do pli
npli = ipli
! Start writing the set of polyline(s).
KMOD = MAX(1,NSH/100)
write(cdigits, '(i1)') int(floor(log10(dble(npli))+1)) ! nr of digits in npli
if (jaNCLAN .or. jaZSH) then
ncol = 3
else
ncol = 2
end if
do ipli=1,npli
if (ipli <= nnam) then
name = names(ipli)
else
name = ' '
end if
! Generate 'L00x' name if empty
if (len_trim(name) <= 0) then
write(name, '(A1,I'//cdigits//'.'//cdigits//')') 'L', ipli
IF (MBNA > 0) THEN
rec = ' '
write(rec , '(A1,I'//cdigits//'.'//cdigits//')') '"', ipli
L = len_trim(rec)
write(rec(L+1:) , '(A)') '","",'
L = len_trim(rec)
write(rec(L+1:) , '(i10)') -(iend(ipli)-istart(ipli)+1)
write(mbna,'(a)') rec
ENDIF
endif
WRITE(MPOL,'(A)') trim(name)
WRITE(MPOL,'(I6,I6)') iend(ipli)-istart(ipli)+1, ncol
! rec = '"324","",-2
DO I = istart(ipli),iend(ipli)
IF (jaNCLAN) THEN
WRITE(MPOL,'(2F15.6,I5)') XSH(I), YSH(I), NCLAN(I)
elseif (jaZSH) then
WRITE(MPOL,'(3F15.6)') XSH(I), YSH(I), ZSH(I)
ELSE
WRITE(MPOL,'(2F15.6)') XSH(I), YSH(I)
IF (MBNA > 0) WRITE(Mbna,'(2F15.6)') XSH(I), YSH(I)
ENDIF
IF (MOD(I,KMOD) .EQ. 0) THEN
CALL READYY(' ',MIN( 1d0,dble(I)/MAX(1,NSH) ) )
ENDIF
END DO ! pts of one polyline
end do ! all polylines
deallocate(istart, iend)
11 CALL READYY(' ',-1d0)
call doclose (MPOL)
IF (MBNA > 0) call doclose (MBNA)
RETURN
END
SUBROUTINE PUTAR (XR,X,MMAX)
implicit none
integer :: i
integer :: mmax
double precision :: x
double precision :: xr
! DE EERSTE IN DE TWEEDE
DIMENSION XR(MMAX) , X(MMAX)
DO 10 I = 1,MMAX
X(I) = XR(I)
10 CONTINUE
RETURN
END
SUBROUTINE DPUTAR (XR,X,MMAX)
implicit none
integer :: i
integer :: mmax
! DE EERSTE IN DE TWEEDE
DOUBLE PRECISION XR(MMAX) , X(MMAX)
DO 10 I = 1,MMAX
X(I) = XR(I)
10 CONTINUE
RETURN
END
SUBROUTINE IPUTAR (IXR,IX,MMAX)
implicit none
integer :: i
integer :: ix
integer :: ixr
integer :: mmax
! DE EERSTE IN DE TWEEDE
DIMENSION IXR(MMAX) , IX(MMAX)
DO 10 I = 1,MMAX
IX(I) = IXR(I)
10 CONTINUE
RETURN
END
SUBROUTINE PUTARR (XR,X,MMAX,NMAX)
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
double precision :: x
double precision :: xr
! DE EERSTE IN DE TWEEDE
DIMENSION XR(MMAX,NMAX), X(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
X(I,J) = XR(I,J)
10 CONTINUE
RETURN
END
SUBROUTINE IPUTARR(XR,X,MMAX,NMAX)
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
! DE EERSTE IN DE TWEEDE
INTEGER XR(MMAX,NMAX), X(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
X(I,J) = XR(I,J)
10 CONTINUE
RETURN
END
SUBROUTINE MINMAX( X, MXLAN, XMIN, XMAX, MAXLAN)
USE M_MISSING
implicit none
integer :: i
integer :: maxlan
integer :: mxlan
double precision :: xmax
double precision :: xmin
double precision :: xx
! BEPAAL MINIMUM EN MAXIMUM VAN EEN EENDIMENSIONALE ARRAY
double precision :: X(MAXLAN)
IF (MXLAN .EQ. 0) THEN
XMIN = 0
XMAX = 0
RETURN
ENDIF
XMIN = 10D20
XMAX = -10D20
DO 10 I = 1,MXLAN
XX = X(I)
IF (XX .NE. dmiss) THEN
XMIN = MIN(XMIN,XX)
XMAX = MAX(XMAX,XX)
ENDIF
10 CONTINUE
IF (XMIN .EQ. 10D20) XMIN = 0
IF (XMAX .EQ.-10D20) XMAX = 0
RETURN
END
SUBROUTINE DMINMAX( X, MXLAN, XMIN, XMAX, MAXLAN)
USE M_MISSING
implicit none
integer :: i
integer :: maxlan
integer :: mxlan
double precision :: xmax
double precision :: xmin
double precision :: xx
! BEPAAL MINIMUM EN MAXIMUM VAN EEN EENDIMENSIONALE ARRAY
DOUBLE PRECISION, intent(inout) :: X(MAXLAN)
IF (MXLAN .EQ. 0) THEN
XMIN = 0
XMAX = 0
RETURN
ENDIF
XMIN = 10D20
XMAX = -10D20
DO 10 I = 1,MXLAN
XX = X(I)
IF (XX .NE. dmiss) THEN
XMIN = MIN(XMIN,XX)
XMAX = MAX(XMAX,XX)
ENDIF
10 CONTINUE
IF (XMIN .EQ. 10D20) XMIN = 0
IF (XMAX .EQ.-10D20) XMAX = 0
RETURN
END
SUBROUTINE DMINMX2( X, XMIN, XMAX, MC,NC,MMAX,NMAX)
USE M_MISSING
implicit none
integer :: i
integer :: j
integer :: mc
integer :: mmax
integer :: nc
integer :: nmax
double precision :: xmax
double precision :: xmin
double precision :: xx
! BEPAAL MINIMUM EN MAXIMUM VAN EEN TWEEDIMENSIONALE ARRAY
DOUBLE PRECISION :: X(MMAX,NMAX)
IF (MC .EQ. 0 .OR. NC .EQ. 0) THEN
XMIN = 0
XMAX = 0
RETURN
ENDIF
XMIN = 10D20
XMAX = -10D20
DO 10 I = 1,MC
DO 10 J = 1,NC
XX = X(I,J)
IF (XX .NE. DXYMIS) THEN
XMIN = MIN(XX,XMIN)
XMAX = MAX(XX,XMAX)
ENDIF
10 CONTINUE
IF (XMIN .EQ. 10D20) XMIN = 0
IF (XMAX .EQ.-10D20) XMAX = 0
RETURN
END
SUBROUTINE MINMXI(IS,MMAX,MINI,MAXI)
USE M_MISSING
implicit none
integer :: i
integer :: maxi
integer :: mini
integer :: mmax
INTEGER IS(MMAX)
MAXI = -999999
MINI = 999999
DO 10 I = 1,MMAX
IF (IS(I) .NE. dmiss) THEN
MAXI = MAX(MAXI,IS(I))
MINI = MIN(MINI,IS(I))
ENDIF
10 CONTINUE
RETURN
END
SUBROUTINE QNREADERROR(W1,W2,MINP)
use unstruc_files
implicit none
integer :: minp
CHARACTER W1*(*),W2*(*)
CALL QNERROR(W1,W2,' IN FILE '//FILENAMES(MINP))
END
SUBROUTINE QNEOFERROR(MINP)
USE unstruc_files
implicit none
integer :: minp
CALL QNERROR('UNEXPECTED END OF FILE IN ',FILENAMES(MINP),' ')
END
SUBROUTINE ININUMBERS()
USE M_MISSING
implicit none
COMMON /NUMBERS/ PI, DG2RD, RD2DG, RA
double precision :: pi, dg2rd, rd2dg, ra
RA = 6370000d0
! RA = dble(6378000.0) DIT IN MEESTE ANDERE LITERATUUR
PI = acos(-1d0)
DG2RD = PI/180d0
RD2DG = 180d0/PI
RETURN
END
LOGICAL FUNCTION INVIEW(X,Y)
USE M_MISSING
use m_wearelt
implicit none
double precision :: x
double precision :: y
! ZIT IK IN ZOOMGEBIED? NULLEN EN DEFAULTS NIET
IF ( X .NE. dmiss .AND. &
X .GT. X1 .AND. X .LT. X2 .AND. &
Y .GT. Y1 .AND. Y .LT. Y2 ) THEN
INVIEW = .TRUE.
ELSE
INVIEW = .FALSE.
ENDIF
RETURN
END
LOGICAL FUNCTION DINVIEW(XD,YD,ZD)
implicit none
double precision :: x
double precision :: y
double precision :: z
LOGICAL INVIEW
DOUBLE PRECISION XD,YD,ZD
CALL DRIETWEE(XD,YD,ZD,X,Y,Z)
DINVIEW = INVIEW(X,Y)
RETURN
END
SUBROUTINE ZEROLAN( KEY)
use M_landboundary
use M_polygon
USE M_MISSING
implicit none
integer :: i
integer :: inhul
integer :: istart
integer :: ja
integer :: k
integer :: key
integer :: mxol
integer :: ntot
KEY = 3
IF (NPL .LE. 2) THEN
CALL CONFRM('NO POLYON, SO DELETE all BOUNDARY POINTS ? ',JA)
IF (JA .EQ. 0) THEN
KEY = 0
RETURN
ENDIF
! CALL SAVESAM()
DO 5 I = 1,MXLAN
XLAN(I) = dmiss
YLAN(I) = dmiss
ZLAN(I) = dmiss
NCLAN(I) = 0
5 CONTINUE
MXLAN = 0
RETURN
ENDIF
! CALL SAVESAM()
INHUL = -1
DO 10 I = 1,MXLAN
CALL DBPINPOL( XLAN(I), YLAN(I), INHUL)
IF (INHUL .EQ. 1) THEN
XLAN(I) = dmiss
YLAN(I) = dmiss
ZLAN(I) = dmiss
NCLAN(I) = 0
ENDIF
10 CONTINUE
K = 0
MXOL = MXLAN
ISTART = 0
NTOT = 0
DO 20 I = 1,MXLAN
IF (XLAN(I) .NE. dmiss) THEN
ISTART = 1
K = K + 1
XLAN(K) = XLAN(I)
YLAN(K) = YLAN(I)
ZLAN(K) = ZLAN(I)
NCLAN(K) = NCLAN(I)
ELSE IF (ISTART .EQ. 1) THEN
K = K + 1
XLAN(K) = dmiss
YLAN(K) = dmiss
ZLAN(K) = dmiss
NCLAN(K) = 0
ISTART = 0
ENDIF
20 CONTINUE
MXLAN = K
DO 30 I = MXLAN+1,MXOL
XLAN(I) = dmiss
YLAN(I) = dmiss
ZLAN(I) = dmiss
NCLAN(I) = 0
30 CONTINUE
RETURN
END
subroutine checkdislin()
use m_polygon
use m_sferic
implicit none
integer :: ja
integer :: jashow
integer :: jmouse
double precision :: xa
double precision :: xlc
double precision :: ya
double precision :: ylc
COMMON /LOCATORA/ XLC,YLC,XA,YA,JMOUSE,JASHOW
double precision :: dis, xn, yn
if (npl >= 2) then
call DLINEDIS(xlc, ylc ,Xpl(1),ypl(1), xpl(2), ypl(2) ,JA,DIS,XN,YN)
endif
call DLINEDIS(1d0,0d0,0d0,0d0,1d0,1d0, JA, DIS,XN,YN)
dis = 0.5d0*sqrt(2d0)
call DLINEDIS(1d0,0.5d0,0d0,0d0,1d0,1d0, JA, DIS,XN,YN)
dis = 0.25d0*sqrt(2d0)
jsferic = 1
call DLINEDIS(4d0,60d0,3d0,60d0,4d0,61d0, JA, DIS,XN,YN)
call DLINEDIS(4d0,60.8d0,3d0,60d0,4d0,61d0, JA, DIS,XN,YN)
call rcirc( xn, yn )
end subroutine checkdislin
subroutine DISPFORM (value,fmt)
implicit none
integer :: n1
integer :: n2
double precision :: value
character fmt*(*)
fmt='(f9.3)'
if (value .eq. 0.0) then
fmt='(f9.5)'
return
endif
n1 = int(log10(abs(value)))
if (n1 .le. 6 .and. n1 .gt. 0) then
n2 = min(9,n1 + 3)
write (fmt(5:5),'(i1)') 9 - n2
else if (n1 .ge. -5 .and. n1 .lt. 0) then
write (fmt(5:5),'(i1)') 6
else if ( n1 .eq. 0) then
write (fmt(5:5),'(i1)') 6
else
fmt ='(e9.3)'
endif
return
end
subroutine DISPFORMscale0(value,fmt)
implicit none
integer :: n1
integer :: n2
integer :: ndec
double precision :: scalesize
double precision :: value
double precision :: xsc
double precision :: ysc
character fmt*(*)
COMMON /SCALEPOS/ XSC,YSC,SCALESIZE,NDEC
fmt='(f9.3)'
if (value .eq. 0.0) then
fmt='(f3.1)'
return
endif
n1 = int(log10(abs(value)))
if (n1 .le. 6 .and. n1 .gt. 0) then
n2 = min(9,n1 + 3)
write (fmt(5:5),'(i1)') 9 - n2
else if (n1 .ge. -5 .and. n1 .lt. 0) then
write (fmt(5:5),'(i1)') 6
else if ( n1 .eq. 0) then
write (fmt(5:5),'(i1)') 6
else
fmt ='(e9.3)'
endif
IF (NDEC .GT. 0) write (fmt(5:5),'(i1)') NDEC
return
end
SUBROUTINE STOPINT()
use unstruc_files
use unstruc_netcdf, only: unc_closeall
use m_partitioninfo
implicit none
CALL ISCREENCLOSE()
call unc_closeall()
call close_all_files()
if ( jampi.eq.1 ) then
! finalize before exit
call partition_finalize()
end if
! SPvdP: close dia-file
if ( mdia.gt.0 .and. mdia.lt.maxnum ) then
close(mdia)
mdia = 0
end if
STOP
END
SUBROUTINE START()
use M_dimens
USE M_DEVICES
use unstruc_files
use unstruc_startup
use unstruc_version_module, only : unstruc_basename
use unstruc_display, only : jaGUI
implicit none
integer :: infofile
integer :: infohardware
integer :: infoopsystem
integer :: ja
integer :: jmouse
integer :: jscreen
integer :: key
integer :: nlevel
integer :: num
integer :: numclargs
integer :: nwhat
COMMON /HELPNOW/ WRDKEY, NLEVEL
COMMON /KERN3D/ INFOFILE,NAMEGRID,NAMEFIELDI,NAMEFIELDO,GRIDAT
COMMON /MESSAGETOSCREEN/ JSCREEN
CHARACTER NAMEGRID*80,NAMEFIELDI*80,NAMEFIELDO*80,GRIDAT*1
CHARACTER WRDKEY*40
!
WRDKEY = 'PROGRAM PURPOSE'
NLEVEL = 1
JSCREEN = 0
INFOFILE = 0
CALL INIDIA(unstruc_basename)
if ( jaGUI.ne.1 ) return
! initialisatiefiles
CALL initProgram()
if ( jaGUI.ne.1 ) return
! SPvdP: disabled mouse-check for mouseless buildserver
! JMOUSE = INFOHARDWARE(13)
! IF (JMOUSE .EQ. 1) THEN
! CALL QNERROR ('NO MOUSE FOUND',' ',' ')
! NLEVEL = 2
! WRDKEY = 'MOUSE INSTALLATION'
! CALL HELP(WRDKEY,NLEVEL)
! CALL STOPINT()
! ENDIF
CALL FIRSTLIN(MDIA)
WRITE(msgbuf,*) 'MAXIMUM NUMBER OF LINKS : ', LMAX ; call msg_flush()
WRITE(msgbuf,*) 'MAXIMUM NUMBER OF NODES : ', KMAX ; call msg_flush()
WRITE(msgbuf,*) 'RESOLUTION GRAPHICS SCREEN : ', NPX, NPY ; call msg_flush()
WRITE(msgbuf,*) 'RESOLUTION TEXT SCREEN : ', IWS, IHS ; call msg_flush()
WRITE(msgbuf,*) 'NUMBER OF COLOURS AVAILABLE : ', NCOLR ; call msg_flush()
15 CONTINUE
NUMCLARGS = INFOOPSYSTEM(2)
IF (NUMCLARGS .GT. 0 .OR. INFOFILE .EQ. 1) RETURN
KEY = 0
JA = 2
CALL MENUH(JA,NUM,NWHAT)
CALL BOTLIN(JA,0,KEY)
IF (KEY .GE. 24 .AND. KEY .LE. 26) THEN
CALL FKEYS(KEY)
GOTO 15
ENDIF
RETURN
END
SUBROUTINE REACLASSES(MINP,CLASS,NUMQ,NUMCLASS)
USE M_MISSING
implicit none
integer :: ja
integer :: k
integer :: l
integer :: minp
integer :: numc
integer :: numclass
integer :: numq
CHARACTER REC*132
double precision :: CLASS(NUMQ,NUMCLASS)
LOGICAL EMPTY
CALL MISARR(CLASS,NUMQ,NUMCLASS)
CALL ZOEKAL(MINp,REC,'TEKAL BLOCK INDICATORS',JA)
if (ja .ne. 1) return
L = 1
10 CONTINUE
READ(MINP,'(A)',END = 999) REC
IF (REC(1:1) .NE. ' ' .OR. EMPTY(REC) ) THEN
NUMC = L - 1
RETURN
ENDIF
READ(REC,*,ERR = 888) (CLASS(K,L),K=1,NUMQ)
L = L + 1
GOTO 10
888 CALL READERROR('READING 5 REALS BUT GETTING',REC,MINP)
999 CALL EOFERROR(MINP)
END
SUBROUTINE READARCINFOHEADER(MINP,MMAX,NMAX,X0,Y0,DX,DY,RMIS)
implicit none
double precision :: dx, dy
integer :: jacornerx
integer :: jacornery
integer :: minp
integer :: mmax
integer :: nmax
double precision :: rmis
double precision :: x0
double precision :: y0
double precision :: DumX, DumY
CHARACTER REC*132
10 CONTINUE
READ(MINP,'(A)',END = 100) REC
IF (REC(1:1) .EQ. '*' .OR. REC(2:2) .EQ. '*') GOTO 10
READ(REC(13:),*,ERR = 101) MMAX
READ(MINP,'(A)',END = 100) REC
READ(REC(13:),*,ERR = 102) NMAX
READ(MINP,'(A)',END = 100) REC
READ(REC(13:),*,ERR = 103) X0
call ilowercase(rec)
JACORNERX = 0
IF (INDEX(REC,'cor') .GE. 1) JACORNERX = 1
READ(MINP,'(A)',END = 100) REC
READ(REC(13:),*,ERR = 104) Y0
call ilowercase(rec)
JACORNERy= 0
IF (INDEX(REC,'cor') .GE. 1) JACORNERy = 1
READ(MINP,'(A)',END = 100) REC
READ(REC(13:),*,ERR = 105) DX
! SPvdP: also try to read a DY on the same line
DY = DX
READ(REC(13:),*,END = 107) DumX, DumY
DY = DumY
107 continue
READ(MINP,'(A)',END = 100) REC
READ(REC(13:),*,ERR = 106) RMIS
IF (JACORNERX .EQ. 1) X0 = X0 + DX/2
IF (JACORNERy .EQ. 1) Y0 = Y0 + DX/2
RETURN
100 CONTINUE
CALL EOFERROR(MINP)
101 CALL READERROR('LOOKING FOR NCOLS (ARC-INFO), BUT GETTING',REC,MINP)
102 CALL READERROR('LOOKING FOR NROWS (ARC-INFO), BUT GETTING',REC,MINP)
103 CALL READERROR('LOOKING FOR XLLCORNER (ARC-INFO), BUT GETTING',REC,MINP)
104 CALL READERROR('LOOKING FOR YLLCORNER (ARCINFO), BUT GETTING',REC,MINP)
105 CALL READERROR('LOOKING FOR CELLSIZE (ARCINFO), BUT GETTING',REC,MINP)
106 CALL READERROR('LOOKING FOR MISSING VALUE (ARCINFO), BUT GETTING',REC,MINP)
END
SUBROUTINE READARCINFOBLOCK(MINP,D,MC,NC,RMIS)
USE M_MISSING
implicit none
integer :: i
integer :: j
integer :: mc
integer :: minp
integer :: nc
double precision :: rmis
double precision :: D(MC,NC)
double precision, dimension(MC) :: dline
CHARACTER TEX*16
DO 10 J = NC,1,-1
READ(MINP,*,ERR=101,END=100) (D(I,J),I = 1,MC)
10 CONTINUE
DO 20 I = 1,MC
DO 20 J = 1,NC
IF (D(I,J) .EQ. RMIS) D(I,J) = dmiss
20 CONTINUE
call doclose (MINP)
RETURN
100 CONTINUE
CALL EOFERROR(MINP)
101 CONTINUE
WRITE(TEX,'(2I8)') I,J
CALL READERROR('ERROR READING ARC-INFO BLOCK IN COLNR, ROWNR :',TEX,MINP)
RETURN
END
!> read Arcinfo data and average it into a smaller array
subroutine ReadLargeArcInfoBlock(MINP, Mfile, Nfile, istart, iend, jstart, jend, Marray, Narray, RMIS, istep, jstep, D)
use m_missing
implicit none
integer, intent(in) :: Mfile, Nfile !< arcinfo dimensions
integer, intent(in) :: istart, iend, jstart, jend !< block to be read in file-index numbering
integer, intent(in) :: Marray, Narray !< sample data array dimensions
integer, intent(inout) :: MINP !< input file unit number
double precision, intent(in) :: RMIS !< missing value
integer, intent(out) :: istep, jstep !< subblock sizes
double precision, dimension(Marray, Narray), intent(inout) :: D !< sample data array
double precision, dimension(:), allocatable :: dline
integer :: iarray, jarray, ifile, jfile, ja3
integer :: isub, jsub
integer, dimension(:), allocatable :: num
double precision :: af, dum
integer :: ierror
CHARACTER TEX*16
ierror = 1
open(6)
! compute subblock sizes
! istep = max(Mfile/Marray,1)
! jstep = max(Nfile/Narray,1)
istep = max((iend-istart+1)/Marray,1)
jstep = max((jend-jstart+1)/Narray,1)
! allocate
allocate(dline(Mfile))
allocate(num(Marray))
D = 0d0
call readyy(' ', -1d0)
call readyy('Reading Arcinfo file (press right mouse button to cancel)', 0d0)
! read last lines outside block
ifile = 1
do jfile=Nfile,jstart+Narray*jstep,-1
read(MINP,*,ERR=101,END=100) ! ( dum, ifile=1,Mfile )
! check for right mouse button
call halt3(ja3)
if ( ja3.eq.3 ) then
! fill remaining array elements with DMISS
D = DMISS
goto 1234
end if
af = dble(jfile-Nfile)/dble(jstart-Nfile)
call readyy('Reading Arcinfo file (press right mouse button to cancel)', af)
WRITE(6,'(1H+"Reading Arcinfo file: ", F7.2, "% done")') 1d2*af
end do
do jarray=Narray,1,-1
! read jstep lines and average
dline = RMIS
num = 0
do jsub=jstep,1,-1
jfile = jstart-1 + jsub + jstep*(jarray-1)
! if ( jfile.gt.Nfile ) cycle
if ( jfile.gt.jend ) cycle
if ( jfile.lt.jstart ) cycle
dline=-1
read(MINP,*,ERR=101,END=100) ( dline(ifile), ifile=1,Mfile )
! sum the sample values in a subcell
do iarray=1,Marray
do isub=1,istep
ifile = istart-1 + isub + istep*(iarray-1)
! if ( ifile.gt.Mfile ) cycle
if ( ifile.gt.iend ) cycle
if ( dline(ifile).ne.RMIS ) then
num(iarray) = num(iarray)+1
D(iarray,jarray) = D(iarray,jarray) + dline(ifile)
end if
end do
end do
end do
af = dble(jfile-Nfile)/dble(jstart-Nfile)
call readyy('Reading Arcinfo file (press right mouse button to cancel)', af)
WRITE(6,'(1H+"Reading Arcinfo file: ", F7.2, "% done")') 1d2*af
! divide by the number of samples in a subcell
do iarray=1,Marray
if (num(iarray).gt.0 ) then
D(iarray,jarray) = D(iarray,jarray) / dble(num(iarray))
else
D(iarray,jarray) = DMISS
end if
end do
! check for right mouse button
call halt3(ja3)
if ( ja3.eq.3 ) then
! fill remaining array elements with DMISS
D(1:Marray,jarray-1:1:-1) = DMISS
goto 1234
end if
end do ! do jarray=Narray,1,-1
call readyy(' ', -1d0)
call doclose (MINP)
ierror = 0
goto 1234
100 CONTINUE
CALL EOFERROR(MINP)
101 CONTINUE
WRITE(TEX,'(2I8)') ifile, jfile
CALL READERROR('ERROR READING ARC-INFO BLOCK IN COLNR, ROWNR :',TEX,MINP)
! error handling
1234 continue
! deallocate
deallocate(dline)
deallocate(num)
close(6)
return
end subroutine ReadLargeArcInfoBlock
SUBROUTINE REAARC(MINP,japrompt)
USE M_ARCINFO
use m_polygon
use m_missing
use m_alloc
implicit none
integer :: ierr
integer :: minp
integer, intent(in) :: japrompt !< prompt for step size (1) or not (0)
integer :: istep, jstep, MCfile, NCfile
integer :: istart, iend, jstart, jend !< block to be read in file-index numbering
double precision :: distep, djstep, dsqrtnumcur
logical :: LdirectReadBlock = .false.
integer, parameter :: MAXSAMSIZE = 1000000
CALL READARCINFOHEADER(MINP,MCa,NCa,X0,Y0,DXa,DYa,RMIS)
IF ( ALLOCATED(D) ) THEN
DEALLOCATE(D)
ENDIF
if ( LdirectReadBlock ) then
ALLOCATE ( D(MCa,NCa),STAT=IERR)
CALL AERR('D(MCa,NCa)',IERR,MCa*NCa)
CALL READARCINFOBLOCK (MINP,D,MCa,NCa,RMIS)
else
istep = 1
jstep = 1
ierr = 1
MCfile = MCa
NCfile = NCa
! do while ( ierr.ne.0 )
if ( NPL.le.0 ) then
istart = 1
iend = MCa
jstart = 1
jend = NCa
else ! use selecting polygon for dimensions of block to be read
istart = max(1+int( (minval(xpl(1:NPL), xpl(1:NPL).ne.DMISS)-X0)/DXa ), 1)
iend = min(1+int( (maxval(xpl(1:NPL), xpl(1:NPL).ne.DMISS)-X0)/DXa ), MCa)
jstart = max(1+int( (minval(ypl(1:NPL), ypl(1:NPL).ne.DMISS)-Y0)/DYa ), 1)
jend = min(1+int( (maxval(ypl(1:NPL), ypl(1:NPL).ne.DMISS)-Y0)/DYa ), NCa)
end if
if ( japrompt.eq.1 ) then
! automatic istep, jstep
dsqrtnumcur = sqrt(dble(iend-istart+1))*sqrt(dble(jend-jstart+1))
distep = dsqrtnumcur/sqrt(dble(MAXSAMSIZE))
distep = dble(int(distep+0.5d0))
djstep = distep
if ( distep.gt.1d0 ) then ! only if necessary
call getreal("istep = ", distep)
call getreal("jstep = ", djstep)
end if
istep = max(int(distep),1)
jstep = max(int(djstep),1)
end if
MCa = (iend-istart+1)/istep
NCa = (jend-jstart+1)/jstep
ALLOCATE ( D(MCa,NCa),STAT=IERR)
! CALL AERR('D(MCa,NCa)',IERR,MCa*NCa)
! check for allocation error
if ( IERR.ne.0 ) then
call qnerror('Sample file too large: increase istep and/or jstep', ' ', ' ')
MCA = 0
NCA = 0
! we cannot deallocate erroneously allocated arrays directly and need to reallocate it correctly first
allocate(D(1,1))
deallocate(D)
goto 1234
end if
! end do ! do while ( ierr.ne.0 )
call ReadLargeArcInfoBlock(MINP, MCfile, NCfile, istart, iend, jstart, jend, MCa, NCa, RMIS, istep, jstep, D)
! modife arcinfo module data
! X0 = X0 + dble(istep-1)*0.5d0*DXa
! Y0 = Y0 + dble(jstep-1)*0.5d0*DYa
X0 = X0 + (istart-1)*Dxa + dble(istep-1)*0.5d0*DXa
Y0 = Y0 + (jstart-1)*Dya + dble(jstep-1)*0.5d0*DYa
DXa = dble(istep)*DXa
DYa = dble(jstep)*DYa
end if ! if ( LdirectReadBlock )
! error handling
1234 continue
RETURN
END
SUBROUTINE WRIARCsam(MARC,DP,MMAX,NMAX,MC,NC,X0,Y0,DX,DY,dmiss)
implicit none
double precision :: dmiss
double precision :: dp
double precision :: dx, dy
integer :: i
integer :: j
integer :: marc
integer :: mc
integer :: mmax
integer :: nc
integer :: nmax
double precision :: x0
double precision :: y0
! DIMENSION DP (MMAX,NMAX)
DIMENSION DP (NMAX,MMAX) ! SPvdP: j-index is fastest running in sample arrays
CALL WRITEARCINFOHEADER(MARC,MC,NC,X0,Y0,DX,DY,dmiss)
! DO 10 J = NC,1,-1 ! SPvdP: j-index is already reverse-order in sample arrays
DO 10 J = 1,NC
! WRITE(MARC,'(5000F10.2)') ( DP(I,J),I = 1,MC)
WRITE(MARC,'(5000F10.2)') ( DP(J,I),I = 1,MC) ! SPvdP: j-index is fastest running in sample arrays
10 CONTINUE
RETURN
END
SUBROUTINE WRIARC(MARC,DP,MMAX,NMAX,MC,NC,X0,Y0,DX,DY,dmiss)
implicit none
double precision :: dmiss
double precision :: dp
double precision :: dx, dy
integer :: i
integer :: j
integer :: marc
integer :: mc
integer :: mmax
integer :: nc
integer :: nmax
double precision :: x0
double precision :: y0
DIMENSION DP (MMAX,NMAX)
CALL WRITEARCINFOHEADER(MARC,MC,NC,X0,Y0,DX,DY,dmiss)
DO 10 J = NC,1,-1
WRITE(MARC,'(5000F10.2)') ( DP(I,J),I = 1,MC)
10 CONTINUE
RETURN
END
SUBROUTINE WRITEARCINFOHEADER(MARC,MC,NC,X0,Y0,DX,DY,dmiss)
implicit none
double precision :: dmiss
double precision :: dx, dy
integer :: marc
integer :: mc
integer :: nc
double precision :: x0
double precision :: y0
WRITE(MARC,'(A,I8)') 'ncols ', MC
WRITE(MARC,'(A,I8)') 'nrows ', NC
WRITE(MARC,'(A,F13.3)') 'xllcorner ', X0 - DX/2
WRITE(MARC,'(A,F13.3)') 'yllcorner ', Y0 - DY/2
WRITE(MARC,'(A,2F13.3)') 'cellsize ', DX, DY
WRITE(MARC,'(A,F13.3)') 'NODATA_value ', dmiss
RETURN
END
SUBROUTINE MISARRR(H,NUMQ,MMAX,NMAX)
USE M_MISSING
implicit none
integer :: i
integer :: j
integer :: k
integer :: mmax
integer :: nmax
integer :: numq
double precision :: H(NUMQ,MMAX,NMAX)
DO 10 K = 1,NUMQ
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
H(K,I,J) = dmiss
10 CONTINUE
RETURN
END
SUBROUTINE MISARR(H,MMAX,NMAX)
USE M_MISSING
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
double precision :: H(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
H(I,J) = dmiss
10 CONTINUE
RETURN
END
SUBROUTINE MISAR(H,MMAX)
USE M_MISSING
implicit none
integer :: i
integer :: mmax
double precision :: H(MMAX)
DO 10 I = 1,MMAX
H(I) = dmiss
10 CONTINUE
RETURN
END
SUBROUTINE IMISARRR(IH,NUMQ,MMAX,NMAX)
USE M_MISSING
implicit none
integer :: i
integer :: j
integer :: k
integer :: mmax
integer :: nmax
integer :: numq
INTEGER IH(NUMQ,MMAX,NMAX)
DO 10 K = 1,NUMQ
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
IH(K,I,J) = dmiss
10 CONTINUE
RETURN
END
SUBROUTINE IMISARR(IH,MMAX,NMAX)
USE M_MISSING
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
INTEGER IH(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
IH(I,J) = dmiss
10 CONTINUE
RETURN
END
SUBROUTINE IMISAR(IH,MMAX)
USE M_MISSING
implicit none
integer :: i
integer :: mmax
INTEGER IH(MMAX)
DO 10 I = 1,MMAX
IH(I) = dmiss
10 CONTINUE
RETURN
END
SUBROUTINE IMISAR2(IH,MMAX)
USE M_MISSING
implicit none
integer :: i
integer :: mmax
INTEGER*2 IH(MMAX)
DO 10 I = 1,MMAX
IH(I) = INT(dmiss)
10 CONTINUE
RETURN
END
SUBROUTINE TMISARR(H,MMAX,NMAX)
USE M_MISSING
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
double precision :: H(NMAX,-1:MMAX+1)
DO 10 I = -1,MMAX+1
DO 10 J = 1,NMAX
H(J,I) = dmiss
10 CONTINUE
RETURN
END
! PLUS QNRGF
SUBROUTINE NUL2(N1,NSMAX)
implicit none
integer :: i
integer :: nsmax
INTEGER*2 N1(NSMAX)
DO 30 I = 1,NSMAX
N1(I) = 0
30 CONTINUE
RETURN
END
SUBROUTINE NULAR (X, MMAX)
implicit none
integer :: i
integer :: mmax
double precision :: x
DIMENSION X(MMAX)
DO 10 I = 1,MMAX
X(I) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE DNULAR (X, MMAX)
implicit none
integer :: i
integer :: mmax
DOUBLE PRECISION X(MMAX)
DO 10 I = 1,MMAX
X(I) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE DNULARR(X, MMAX, NMAX)
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
DOUBLE PRECISION X(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
X(I,J) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE NULARR(X, MMAX, NMAX)
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
double precision :: x
DIMENSION X(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
X(I,J) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE NULARRR(X, MMAX, NMAX, LMAX)
implicit none
integer :: i
integer :: j
integer :: l
integer :: lmax
integer :: mmax
integer :: nmax
double precision :: x
DIMENSION X(MMAX,NMAX,LMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
DO 10 L = 1,LMAX
X(I,J,L) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE INULAR (X, MMAX)
implicit none
integer :: i
integer :: mmax
INTEGER X(MMAX)
DO 10 I = 1,MMAX
X(I) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE INULARR(X, MMAX, NMAX)
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
INTEGER X(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
X(I,J) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE INULARRR(X, MMAX, NMAX, LMAX)
implicit none
integer :: i
integer :: j
integer :: l
integer :: lmax
integer :: mmax
integer :: nmax
INTEGER X(MMAX,NMAX,LMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
DO 10 L = 1,LMAX
X(I,J,L) = 0
10 CONTINUE
RETURN
END
SUBROUTINE INULAR2(X, MMAX)
implicit none
integer :: i
integer :: mmax
INTEGER*2 X(MMAX)
DO 10 I = 1,MMAX
X(I) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE INULARR2(X, MMAX, NMAX)
implicit none
integer :: i
integer :: j
integer :: mmax
integer :: nmax
INTEGER*2 X(MMAX,NMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
X(I,J) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE INULARRR2(X, MMAX, NMAX, LMAX)
implicit none
integer :: i
integer :: j
integer :: l
integer :: lmax
integer :: mmax
integer :: nmax
INTEGER*2 X(MMAX,NMAX,LMAX)
DO 10 I = 1,MMAX
DO 10 J = 1,NMAX
DO 10 L = 1,LMAX
X(I,J,L) = 0d0
10 CONTINUE
RETURN
END
SUBROUTINE RE0RCINFODIMENSIONS(MINP,MMAX,NMAX,DX,X0,Y0)
implicit none
double precision :: dx
integer :: ja
integer :: larc
integer :: marc
integer :: mc
integer :: minp
integer :: mmax
integer :: nc
integer :: nmax
double precision :: rmis
double precision :: x0
double precision :: y0
CHARACTER REC*132, FILENAME*80
REWIND(MINP)
CALL ZOEKJA(MINP,REC,'ARC-INFO',JA)
IF (JA .EQ. 1) THEN
LARC = INDEX(REC,'ARC')
READ(REC(LARC+9:),'(A)',ERR = 888) FILENAME
CALL OLDFIL(MARC,FILENAME)
CALL READARCINFOHEADER(MARC,MC,NC,X0,Y0,DX,dx,RMIS)
call doclose(MARC)
MMAX = MC
NMAX = NC
CALL MESSAGE('MAIN DIMENSIONS PLUS GRID PARAMETERS HAVE', &
'BEEN READ FROM ARC-INFO FILE:', FILENAME)
ELSE
CALL ERROR('NEITHER MAIN DIMENSIONS NOR ARC-INFO FILE FOUND,', &
'SPECIFY MAIN DIMENSIONS BY KEYWORD:', &
'MAIN DIMENSIONS OR BY ARC-INFO FILE')
ENDIF
RETURN
888 CALL ERROR('LOOKING FOR ARC-INFO FILENAME, BUT GETTING:',REC,' ')
END
subroutine dateandtimenow(iyear, month, iday, ihour, minute, isecnd)
implicit none
integer, intent(out) :: iyear, month, iday, ihour, minute, isecnd
! integer, intent(out), optional :: imsec
! character(len=5), intent(out), optional :: zone
character(len=8 ) :: dat
character(len=10) :: tim
character(len=5) :: zone
integer :: imsec
integer :: values(8)
call date_and_time(dat,tim,zone,values)
iyear = values(1)
month = values(2)
iday = values(3)
ihour = values(5)
minute = values(6)
isecnd = values(7)
imsec = values(8)
end subroutine dateandtimenow
SUBROUTINE DATUM(DATE)
implicit none
integer :: iyear, month, iday, ihour, minute, isecnd
CHARACTER DATE*20
! 1 4 7 11 14 17
DATE = 'hh:mm:ss, dd-mm-yyyy'
call dateandtimenow(iyear, month, iday, ihour, minute, isecnd)
WRITE(DATE( 1:2 ),'(I2.2)') IHOUR
WRITE(DATE( 4:5 ),'(I2.2)') MINUTE
WRITE(DATE( 7:8 ),'(I2.2)') ISECND
WRITE(DATE(11:12),'(I2.2)') IDAY
WRITE(DATE(14:15),'(I2.2)') MONTH
WRITE(DATE(17:20),'(I4)') IYEAR
RETURN
END
SUBROUTINE DATUM2(DATE)
use unstruc_display, only : jadatetime
implicit none
integer :: iyear, month, iday, ihour, minute, isecnd
CHARACTER DATE*20
character(len=1), external :: get_dirsep
! 1 4 7 11 14 17
if (jadatetime == 0) then
DATE = get_dirsep()
else
DATE = '_yymmddhhmmss'//get_dirsep()
call dateandtimenow(iyear, month, iday, ihour, minute, isecnd)
WRITE(DATE(2:3),'(I2.2)') IYEAR - 2000
WRITE(DATE(4:5),'(I2.2)') month
WRITE(DATE(6:7),'(I2.2)') iday
WRITE(DATE(8:9),'(I2.2)') Ihour
WRITE(DATE(10:11),'(I2.2)') minute
WRITE(DATE(12:13),'(I2.2)') isecnd
endif
RETURN
END
integer FUNCTION USECP()
implicit none
integer :: ihour
integer :: isecnd
integer :: minute
integer :: iy,im,id
call dateandtimenow(iy,im,id,IHOUR,MINUTE,ISECND)
USECP = 3600*IHOUR+60*MINUTE+ISECND
RETURN
END
SUBROUTINE READXYMIS(MINP)
USE M_MISSING
implicit none
integer :: ja
integer :: l
integer :: minp
! snelheidsdrempel
CHARACTER REC*132, TEX*8
CALL ZOEKAL(MINP,REC,'MISSING VALUE XY',JA)
XYMIS = 0d0
IF (JA .EQ. 1) THEN
L = INDEX(REC,'XY') + 4
READ(REC(L:),*,ERR = 888) XYMIS
WRITE(TEX,'(F8.3)') XYMIS
CALL MESSAGE('MISSING VALUE XY = ', TEX,' ')
ENDIF
RETURN
888 CALL READERROR('READING MISSING VALUE XY, BUT GETTING', REC, MINP)
END
SUBROUTINE INDEKSI(N,NARRIN,INDX)
implicit none
integer :: i
integer :: indx
integer :: indxt
integer :: ir
integer :: j
integer :: l
integer :: n
integer :: narrin
double precision :: q
DIMENSION NARRIN(N),INDX(N)
DO 11 J=1,N
INDX(J)=J
11 CONTINUE
L=N/2+1
IR=N
10 CONTINUE
IF(L.GT.1)THEN
L=L-1
INDXT=INDX(L)
Q=NARRIN(INDXT)
ELSE
INDXT=INDX(IR)
Q=NARRIN(INDXT)
INDX(IR)=INDX(1)
IR=IR-1
IF(IR.EQ.1)THEN
INDX(1)=INDXT
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF(J.LE.IR)THEN
IF(J.LT.IR)THEN
IF(NARRIN(INDX(J)).LT.NARRIN(INDX(J+1)))J=J+1
ENDIF
IF(Q.LT.NARRIN(INDX(J)))THEN
INDX(I)=INDX(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
INDX(I)=INDXT
GO TO 10
END
!
SUBROUTINE INDEKS(N,ARRIN,INDX)
implicit none
double precision :: arrin
integer :: i
integer :: indx
integer :: indxt
integer :: ir
integer :: j
integer :: l
integer :: n
double precision :: q
DIMENSION ARRIN(N),INDX(N)
DO 11 J=1,N
INDX(J)=J
11 CONTINUE
L=N/2+1
IR=N
10 CONTINUE
IF(L.GT.1)THEN
L=L-1
INDXT=INDX(L)
Q=ARRIN(INDXT)
ELSE
INDXT=INDX(IR)
Q=ARRIN(INDXT)
INDX(IR)=INDX(1)
IR=IR-1
IF(IR.EQ.1)THEN
INDX(1)=INDXT
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF(J.LE.IR)THEN
IF(J.LT.IR)THEN
IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
ENDIF
IF(Q.LT.ARRIN(INDX(J)))THEN
INDX(I)=INDX(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
INDX(I)=INDXT
GO TO 10
END
SUBROUTINE BILIN(X, Y, Z, XP, YP, ZP)
implicit none
double precision :: r1
double precision :: r2
double precision :: x1
double precision :: x2
double precision :: xa
double precision :: xp
double precision :: xr
double precision :: xrm
double precision :: y1
double precision :: y2
double precision :: ya
double precision :: yp
double precision :: yr
double precision :: yrm
double precision :: zp
double precision :: X(4), Y(4), Z(4)
! Bepaal relatieve ligging in cel.
! Twee coordinaten (xr,yr) van 0 tot 1.
! Bilineaire interpolatie
X1 = X(2) - X(1)
Y1 = Y(2) - Y(1)
R1 = X1*X1 + Y1*Y1
X2 = X(4) - X(1)
Y2 = Y(4) - Y(1)
R2 = X2*X2 + Y2*Y2
XA = XP - X(1)
YA = YP - Y(1)
XR = (X1*XA + Y1*YA)/R1
YR = (X2*XA + Y2*YA)/R2
XRM = 1 - XR
YRM = 1 - YR
ZP = XRM*YRM*Z(1) + XR*YRM*Z(2) + XR*YR*Z(3) + XRM*YR*Z(4)
RETURN
END
SUBROUTINE GETDIMENSIONS(MXD,NXD,MXLN,NSX)
implicit none
integer :: mout
integer :: mxd
integer :: mxln
integer :: nsx
integer :: nxd
CHARACTER GETAL*100
LOGICAL THISISANUMBER, JAWEL
MXD = 500 ! ROOSTERS EN SPLINES M-RICHTING
NXD = 500 ! ROOSTERS EN SPLINES N-RICHTING
MXLN = 100000 ! land boundary
NSX = 100000 ! SAMPLES
GETAL = ' '
CALL get_command_argument(1,GETAL)
IF (THISISANUMBER(GETAL)) READ(GETAL,*) MXD
GETAL = ' '
CALL get_command_argument(2,GETAL)
IF (THISISANUMBER(GETAL)) READ(GETAL,*) NXD
GETAL = ' '
CALL get_command_argument(3,GETAL)
IF (THISISANUMBER(GETAL)) READ(GETAL,*) MXLN
GETAL = ' '
CALL get_command_argument(4,GETAL)
IF (THISISANUMBER(GETAL)) READ(GETAL,*) NSX
INQUIRE (FILE = 'rgfdim', EXIST = JAWEL)
IF (JAWEL) THEN
MOUT = 10
call oldfil(MOUT,'rgfdim')
READ (MOUT,*,ERR=999) MXD,NXD,MXLN,NSX
call doclose(MOUT)
ENDIF
999 CONTINUE
RETURN
END
SUBROUTINE BILINXY(X, Y, XZ, YZ, XP, YP, XP2, YP2, INI)
implicit none
double precision :: c
integer :: i
integer :: ini
integer :: japarallel
double precision, SAVE :: A(4,4), BX(4), BY(4)
double precision :: X(4), Y(4), XZ(4), YZ(4), XP, YP, XP2, YP2
INTEGER, SAVE :: INX(4)
! (Zi = AXi + BYi + CXiYi + Di ,i=1,4)
! Coefficienten in A, rechterlid in B, opl met LU-decompositie
IF (INI .EQ. 1) THEN
DO I = 1,4
A(I,1) = X(I)-X(1)
A(I,2) = Y(I)-Y(1)
A(I,3) = (Y(I)-Y(1))*(X(I)-X(1))
A(I,4) = 1
BX(I) = XZ(I)
BY(I) = YZ(I)
ENDDO
CALL LUDCMP(A,4,4,INX,C,JAPARALLEL)
IF (JAPARALLEL .EQ. 1) THEN
CALL qnerror('Problem in Ludcmp',' ',' ')
INI = -1
RETURN
ENDIF
CALL LUBKSB(A,4,4,INX,BX)
CALL LUBKSB(A,4,4,INX,BY)
ENDIF
XP2 = (XP-X(1))*BX(1) + (YP-Y(1))*BX(2) + (XP-X(1))*(YP-Y(1))*BX(3) + BX(4)
YP2 = (XP-X(1))*BY(1) + (YP-Y(1))*BY(2) + (XP-X(1))*(YP-Y(1))*BY(3) + BY(4)
RETURN
END
SUBROUTINE LUDCMP(A,N,NP,INDX,D,JAPARALLEL)
implicit none
double precision :: a
double precision :: aamax
double precision :: d
double precision :: dum
integer :: i
integer :: imax
integer :: indx
integer :: j
integer :: japarallel
integer :: k
integer :: n
integer :: np
integer :: nx
double precision :: sum
double precision :: tiny
double precision :: vv
PARAMETER (NX=4,TINY=1d-20)
DIMENSION A(NP,NP),INDX(N),VV(NX)
JAPARALLEL = 0
D=1.
DO 12 I=1,N
AAMAX=0.
DO 11 J=1,N
IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11 CONTINUE
IF (AAMAX .EQ. 0) THEN
JAPARALLEL = 1
RETURN
ENDIF
VV(I)=1./AAMAX
12 CONTINUE
DO 19 J=1,N
IF (J.GT.1) THEN
DO 14 I=1,J-1
SUM=A(I,J)
IF (I.GT.1)THEN
DO 13 K=1,I-1
SUM=SUM-A(I,K)*A(K,J)
13 CONTINUE
A(I,J)=SUM
ENDIF
14 CONTINUE
ENDIF
AAMAX=0.
DO 16 I=J,N
SUM=A(I,J)
IF (J.GT.1)THEN
DO 15 K=1,J-1
SUM=SUM-A(I,K)*A(K,J)
15 CONTINUE
A(I,J)=SUM
ENDIF
DUM=VV(I)*ABS(SUM)
IF (DUM.GE.AAMAX) THEN
IMAX=I
AAMAX=DUM
ENDIF
16 CONTINUE
IF (J.NE.IMAX)THEN
DO 17 K=1,N
DUM=A(IMAX,K)
A(IMAX,K)=A(J,K)
A(J,K)=DUM
17 CONTINUE
D=-D
VV(IMAX)=VV(J)
ENDIF
INDX(J)=IMAX
IF(J.NE.N)THEN
IF(A(J,J).EQ.0d0)A(J,J)=TINY
DUM=1./A(J,J)
DO 18 I=J+1,N
A(I,J)=A(I,J)*DUM
18 CONTINUE
ENDIF
19 CONTINUE
IF(A(N,N).EQ.0d0)A(N,N)=TINY
RETURN
END
SUBROUTINE LUBKSB(A,N,NP,INDX,B)
implicit none
double precision :: a
double precision :: b
integer :: i
integer :: ii
integer :: indx
integer :: j
integer :: ll
integer :: n
integer :: np
double precision :: sum
DIMENSION A(NP,NP),INDX(N),B(N)
II=0
DO 12 I=1,N
LL=INDX(I)
SUM=B(LL)
B(LL)=B(I)
IF (II.NE.0)THEN
DO 11 J=II,I-1
SUM=SUM-A(I,J)*B(J)
11 CONTINUE
ELSE IF (SUM.NE.0d0) THEN
II=I
ENDIF
B(I)=SUM
12 CONTINUE
DO 14 I=N,1,-1
SUM=B(I)
IF(I.LT.N)THEN
DO 13 J=I+1,N
SUM=SUM-A(I,J)*B(J)
13 CONTINUE
ENDIF
B(I)=SUM/A(I,I)
14 CONTINUE
RETURN
END
SUBROUTINE SETCLASSFILE()
implicit none
integer :: minp
CHARACTER FILNAM*76
FILNAM = '*.cls'
MINP = 0
CALL FILEMENU(MINP,FILNAM)
IF (MINP .GT. 0) THEN
CALL REACLS(MINP)
ENDIF
RETURN
END
SUBROUTINE SYSORLOCALFIL(LUNID,FILNAM,MUSTBE)
use unstruc_messages
use unstruc_files
implicit none
integer :: ifirstchar
integer :: istart
integer :: k1
integer :: k2
integer :: lunid
integer :: mustbe
integer :: istat
CHARACTER FILNAM*76,FULNAM*180
LOGICAL JA
LUNID = 0
INQUIRE(FILE = FILNAM,EXIST = JA)
IF (JA) THEN
CALL OLDFIL(LUNID,FILNAM)
WRITE(msgbuf,'(2A)') 'Using Local File', FILNAM; call msg_flush()
ELSE
FULNAM = PATHDI
ISTART = len_trim(PATHDI) + 1
WRITE(FULNAM(ISTART:),'(A)') FILNAM
K1 = IFIRSTCHAR (FULNAM)
K2 = len_trim(FULNAM)
INQUIRE(FILE= FULNAM(K1:K2),EXIST=JA)
IF (JA) THEN
CALL OLDFIL(LUNID,FULNAM)
call mess(LEVEL_INFO, 'Using Program File ', FULNAM(K1:K2))
ELSE IF (MUSTBE .EQ. 1) THEN
call mess(LEVEL_ERROR, 'Program File ',FULNAM(K1:K2),' Not Found')
ENDIF
ENDIF
RETURN
END
SUBROUTINE REACLS (MLAN) ! DV
implicit none
double precision :: dv
integer :: i
integer :: jaauto
integer :: mlan
integer :: ncols
integer :: nie
integer :: nis
integer :: nrow
integer :: nv
double precision :: val
double precision :: vmax
double precision :: vmin
double precision :: x
! ------------------------------------------------------------------
! LEZEN FILE MET CLASSES
! ------------------------------------------------------------------
COMMON /DEPMAX/ VMAX,VMIN,DV,VAL(256),NCOLS(256),NV,NIS,NIE,JAAUTO
CHARACTER MATR*4
CALL READYY ('READING CLS-FILE',0d0)
! ------------------------------------------------------------------
! EERST LEZEN ALS DUMMY WAARDEN OM TE TESTEN OF ER FOUTEN OPTREDEN
! ------------------------------------------------------------------
10 READ(MLAN,'(A)',END=999) MATR
IF (MATR(1:1) .EQ. '*') GOTO 10
READ (MLAN,*,ERR=999) NROW
NROW = MIN(NROW,30)
DO 20 I = 1,NROW
READ (MLAN,*,ERR=999,END=999) X
20 CONTINUE
! -------------------------------------
! NO ERRORS, SO READ THE VALUES
! -------------------------------------
REWIND (MLAN)
110 READ(MLAN,'(A)',END=999) MATR
IF (MATR(1:1) .EQ. '*') GOTO 110
READ (MLAN,*,ERR=999) NROW
NROW = MIN(NROW,30)
DO 120 I = 1,NROW
READ (MLAN,*) VAL(I)
VMIN = MIN(VMIN,VAL(I))
VMAX = MAX(VMAX,VAL(I))
120 CONTINUE
JAAUTO = 2
NV = NROW
999 CONTINUE
CALL READYY(' ', 1d0)
CALL READYY(' ',-1d0)
call doclose (MLAN)
RETURN
END
subroutine triinterp2(XZ,YZ,BL,NDX,JDLA)
use m_samples
use m_polygon
USE M_MISSING
implicit none
integer :: jdla
integer :: jakdtree
integer :: k
integer :: ndx
DOUBLE PRECISION :: XZ(ndx), YZ(ndx), BL(ndx)
if (ndx < 1) return
jakdtree = 1
if ( MXSAM.gt.0 .and. MYSAM.gt. 0 ) then ! bi-linear interpolation
call bilin_interp(NDX, XZ, YZ, BL)
else ! Delauny
call TRIINTfast(XS,YS,ZS,NS,1,XZ,YZ,BL,Ndx,Xpl,Ypl,Npl,JDLA,jakdtree)
end if
end subroutine triinterp2
SUBROUTINE TRIINTfast(XS,YS,ZS,NS,NDIM,X,Y,Z,NXY,XH,YH,NH,JDLA,jakdtree)
USE M_TRIANGLE
USE M_MISSING
use unstruc_colors
use m_interpolationsettings
! use m_kdtree2, only: ITREE_EMPTY, ITREE_DIRTY
use m_sferic, only: jsferic
use m_kdtree2
use m_flowexternalforcings , only: transformcoef
implicit none
integer, intent(inout) :: jakdtree !< use kdtree (1) or not (0)
double precision :: af
integer :: i, IERR, k, inhul, intri, jdla, jslo
integer :: n, in, nxx, nxy
integer :: NDIM !< sample vector dimension
integer :: naf
integer :: nbf
integer :: ncf
integer :: ncol
integer :: ndraw
integer :: nh
integer :: nrfind
integer :: ns, n2
integer :: idim
double precision :: rd
double precision :: slo(NDIM)
double precision :: xmaxs
double precision :: xmins
double precision :: xp, yp, zp(NDIM), xpmin, xpmax,ypmin, ypmax
double precision :: ymaxs
double precision :: ymins
double precision :: XS(ns), YS(ns), ZS(NDIM,ns), X(nxy), Y(nxy), Z(NDIM,nxy)
double precision :: XH(nh), YH(nh), XL(3),YL(3)
double precision, allocatable :: xx(:), yy(:) , zz(:,:)
integer , allocatable :: ind(:)
double precision, dimension(:), allocatable :: xs1, ys1 ! for store/restore of xs, ys
double precision, dimension(:,:), allocatable :: zs1 ! for store/restore of zs
double precision :: xsmin, ysmin, xsmax, ysmax ! bounding box corner coordinates
integer :: indf(3)
double precision :: wf(3)
integer :: NS1 ! for store/restore of NS
integer :: jadum, ierror
logical :: Ldeleteddata ! for store/restore of xs, ys and zs
type(kdtree_instance) :: sampletree
COMMON /DRAWTHIS/ NDRAW(40)
integer :: numsam, ii, k1, jakdtree2 ! Nabi
double precision :: R2Search, dist2, cof1, cof2
jakdtree2 = jakdtree
if ( jakdtree.eq.1 ) then
! enforce generation of kdtree
treeglob%itreestat = ITREE_EMPTY
end if
! JSLO=1, OOK SLOPES RD4
IF (NS .LE. 2) RETURN
JSLO = 0
Ldeleteddata = .false.
jdla = 1 ! later, maybe remove this if you want to avoid dlauny every time. This is not working now.
IF (JDLA == 1) THEN
IF (Nh > 2) THEN ! polygon size reduction
xpmin = minval(xh(1:nh)) ; xpmax = maxval(xh(1:nh))
ypmin = minval(yh(1:nh)) ; ypmax = maxval(yh(1:nh))
rd = 0.2*(xpmax-xpmin) ; xpmin = xpmin - rd ; xpmax = xpmax + rd
rd = 0.2*(ypmax-ypmin) ; ypmin = ypmin - rd ; ypmax = ypmax + rd
! CALL SAVESAM() ! xs, ys, zs not necessarilly sample data
! store original data
allocate(xs1(NS), ys1(NS), zs1(NDIM,NS))
NS1 = NS
do i=1,NS
xs1(i) = xs(i)
ys1(i) = ys(i)
do k=1,NDIM
zs1(k,i) = zs(k,i)
end do
end do
Ldeleteddata = .true.
n2 = 0
idim = 1 ! DMISS check on first sample dimension only
in = -1
DO K = 1,Ns
if (jins == 1) then
if ( xs(k) > xpmin .and. xs(k) < xpmax .and. ys(k) > ypmin .and. ys(k) < ypmax .and. zs(idim,k).ne.DMISS ) then
n2 = n2 + 1; xs(n2) = xs(k) ; ys(n2) = ys(k) ; zs(1:NDIM,n2) = zs(1:NDIM,k)
endif
else
if (zs(idim,k).ne.DMISS ) then
n2 = n2 + 1; xs(n2) = xs(k) ; ys(n2) = ys(k) ; zs(1:NDIM,n2) = zs(1:NDIM,k)
endif
endif
ENDDO
nxx = 0 ! count net/grid size
DO K = 1,Nxy
if (jins == 1) then
if ( x(k) > xpmin .and. x(k) < xpmax .and. y(k) > ypmin .and. y(k) < ypmax ) then
nxx = nxx + 1
endif
else
nxx = nxx + 1
endif
ENDDO
allocate (xx (nxx), yy(nxx), zz(NDIM,nxx) )
allocate (ind(nxx)) ; ind = 0
nxx = 0 ; in = -1 ! net/grid size reduction
DO K = 1,Nxy
call dbpinpol(x(k), y(k), in)
if (in == 1) then
nxx = nxx + 1
ind(nxx) = k
xx (nxx) = x(k); yy(nxx) = y(k); zz(1:NDIM,nxx) = z(1:NDIM,k)
endif
ENDDO
ELSE
NXX = NXY; N2 = NS
ENDIF
! determine sample bounding box
xsmin = minval(xs,mask=xs.ne.DMISS)
xsmax = maxval(xs,mask=xs.ne.DMISS)
ysmin = minval(ys,mask=ys.ne.DMISS)
ysmax = maxval(ys,mask=ys.ne.DMISS)
ierr = 1
if ( jakdtree.ne.1 ) then
CALL DLAUN(XS,YS,N2,1,IERR)
else
call dlaun(xs,ys,N2,3,ierr) ! generate edgeindex and triedge
end if
if ( IERR.ne.0 ) then
goto 1234
end if
JDLA = 0
ENDIF
IF (NUMTRI .LT. 1) THEN
RETURN
ENDIF
!if ( numtri.lt.40 ) then
! jakdtree=0
!end if
NCOL = 14
IF (Jtekinterpolationprocess .NE. 0) THEN
DO I = 1,NUMTRI
NAF = INDX(1,I)
NBF = INDX(2,I)
NCF = INDX(3,I)
XL(1) = XS(NAF)
XL(2) = XS(NBF)
XL(3) = XS(NCF)
YL(1) = YS(NAF)
YL(2) = YS(NBF)
YL(3) = YS(NCF)
CALL TEKTRI(XL,YL,NCOLDN)
ENDDO
ENDIF
R2search = 0d0
if( transformcoef(6) /= dmiss ) R2search = transformcoef(6)**2
if ( jakdtree2 == 1 .and. R2search.gt.0d0 ) then ! build kd-tree for sample points
sampletree%itreestat = ITREE_EMPTY
call build_kdtree(sampletree, n2,xs,ys,ierror)
if ( ierror /= 0 ) then
if ( sampletree%itreestat /= ITREE_EMPTY ) call delete_kdtree2(sampletree)
jakdtree2 = 0
end if
end if
KMOD = MAX(1,NXX/100)
DO N = 1,NXx
IF (Jtekinterpolationprocess == 0 .AND. MOD(N,KMOD) == 0 ) THEN
AF = dble(N-1) / dble(NXX)
CALL READYY('TRIANGLE INTERPOLATION',AF)
ENDIF
idim = 1 ! DMISS check on first sample dimension only
if (nh > 2) then
XP = xx(N)
YP = yy(N)
RD = zz(idim,N)
else
XP = x(N)
YP = y(N)
RD = z(idim,N)
endif
INTRI = 0
! For triangulation, only consider points that are inside the sample bounding box
if ( xp.ge.xsmin .and. xp.le.xsmax .and. yp.ge.ysmin .and. yp.le.ysmax ) then
IF (RD .EQ. dmiss) then
if ( jakdtree.eq.1 ) then
jadum = 0
!if ( N.eq.974271 ) jadum = 1
!if ( N.eq.21 ) jadum = 1
!if ( N.eq.263 ) jadum = 1
!if ( N.eq.5839 .or. N.eq.5935 ) jadum = 1
!if ( N.eq.31 .or. N.eq.169 .or. N.eq.360 ) jadum = 1
!if ( N.eq.8081 ) jadum = 1
!if ( N.eq.7797 ) jadum = 1
!if ( N.eq.341 ) jadum = 1
call findtri_kdtree(XP,YP,ZP,XS,YS,ZS,N2,NDIM, &
NRFIND,INTRI,JSLO,SLO,Jtekinterpolationprocess,jadum,ierror,indf,wf)
if ( ierror.ne.0 ) then
! deallocate
if ( treeglob%itreestat.ne.ITREE_EMPTY ) call delete_kdtree2(treeglob)
if ( allocated(imask) ) deallocate(imask)
if ( allocated(LNtri) ) deallocate(LNtri)
! disable kdtree
jakdtree = 0
end if
end if
if ( jakdtree.ne.1 ) then
CALL FINDTRI(XP,YP,ZP,XS,YS,ZS,N2,NDIM, &
NRFIND,INTRI,JSLO,SLO,Jtekinterpolationprocess,indf,wf)
end if
IF (INTRI .EQ. 1) THEN
if (nh > 2) then
Zz(:,N) = ZP
else
z(:,n) = zp
endif
if (jagetwf == 1) then
indxx(1:3,n) = indf(1:3)
wfxx (1:3,n) = wf (1:3)
endif
ENDIF ! (INTRI .EQ. 1)
ENDIF ! (XP .NE. XYMIS .AND. ... .AND. YP .LE. YMAXS )
endif
!!!!!!!!!! give it another try with nearest neighbour opr inverse distance.
if (intri == 0 .and. R2search.gt.0d0) then
if (RD == dmiss) then ! Nabi
if( jakdtree2 == 1 ) then
call make_queryvector_kdtree(sampletree, xp, yp )
numsam = kdtree2_r_count( sampletree%tree, sampletree%qv, R2search )
if ( numsam.gt.0 ) then
call realloc_results_kdtree(sampletree, numsam)
call kdtree2_n_nearest(sampletree%tree,sampletree%qv,numsam,sampletree%results)
end if
do i = 1,idim
ii = 0
cof1 = 0d0
cof2 = 0d0
do k = 1,numsam
k1 = sampletree%results(k)%idx
if( abs( xp - xs(k1) ) < 1d-6 .and. abs( yp - ys(k1) ) < 1d-6 ) then
ii = 1
z(i,n) = zs(i,k1)
exit
endif
dist2 = ( xp - xs(k1) )**2 + ( yp - ys(k1) )**2
cof1 = cof1 + zs(i,k1) / dist2
cof2 = cof2 + 1d0 / dist2
enddo
if( ii == 0 .and. numsam > 0 ) then
z(i,n) = cof1 / cof2
end if
enddo
else
do i = 1,idim
ii = 0
cof1 = 0d0
cof2 = 0d0
do k = 1,n2
if( abs( xp - xs(k) ) < 1d-6 .and. abs( yp - ys(k) ) < 1d-6 ) then
ii = 1
z(i,n) = zs(i,k)
exit
endif
dist2 = ( xp - xs(k) )**2 + ( yp - ys(k) )**2
cof1 = cof1 + zs(i,k) / dist2
cof2 = cof2 + 1d0 / dist2
enddo
if( ii == 0 ) then
z(i,n) = cof1 / cof2
end if
enddo
endif
endif ! RDMISS
end if ! intri ==0
ENDDO ! DO N = 1,NXx
IF (Jtekinterpolationprocess == 0) CALL READYY('TRIANGLE INTERPOLATION',-1d0)
if (nh > 2) then
do k = 1,nxx
do idim=1,NDIM
z(idim,ind(k)) = zz(idim,k)
end do
enddo
deallocate (xx, yy, zz, ind)
! CALL RESTORESAM()
endif
1234 continue
! deallocate
if ( jakdtree.eq.1 ) then
if ( treeglob%itreestat .ne. ITREE_EMPTY ) call delete_kdtree2(treeglob )
if ( sampletree%itreestat .ne. ITREE_EMPTY ) call delete_kdtree2(sampletree)
if ( allocated(imask) ) deallocate(imask)
if ( allocated(LNtri) ) deallocate(LNtri)
end if
! restore original data
if ( Ldeleteddata ) then
NS = NS1
do i=1,NS
xs(i) = xs1(i)
ys(i) = ys1(i)
do k=1,NDIM
zs(k,i) = zs1(k,i)
end do
end do
deallocate(xs1, ys1, zs1)
end if
RETURN
END subroutine triintfast
SUBROUTINE FINDTRI(XP,YP,ZP,XS,YS,ZS,NS,NDIM,NRFIND,INTRI,JSLO,SLO,JATEK, ind, wf)
USE M_MISSING
USE M_TRIANGLE
implicit none
integer :: intri, jatek, jslo
integer :: k, k1, k2, numit
integer :: nrfind, nroldfind, interval
integer :: ns
integer :: NDIM !< sample vector dimension
integer :: idim, ind(3)
double precision :: slo(NDIM)
double precision :: xp
double precision :: xtmax
double precision :: xtmin
double precision :: yp
double precision :: ytmax
double precision :: ytmin
double precision :: zp(NDIM)
double precision :: XS(NS), YS(NS), ZS(NDIM,NS), XT(3),YT(3),ZT(NDIM,3), wf(3)
integer :: ik1, ik2, ik3
DATA NROLDFIND /0/
ZP = dmiss
INTRI = 0
interval = 2 ; numit = 0; nrfind = 0
5 CONTINUE
numit = numit + 1
interval = 5*interval
K1 = MAX(1 ,NROLDFIND-interval)
K2 = MIN(NUMTRI,NROLDFIND+interval)
DO K = K1,K2
ik1 = INDX(1,K)
ik2 = INDX(2,K)
ik3 = INDX(3,K)
XT(1) = XS(ik1)
XT(2) = XS(ik2)
XT(3) = XS(ik3)
YT(1) = YS(ik1)
YT(2) = YS(ik2)
YT(3) = YS(ik3)
XTMAX = MAX(XT(1),MAX( XT(2),XT(3) ) )
YTMAX = MAX(YT(1),MAX( YT(2),YT(3) ) )
XTMIN = MIN(XT(1),MIN( XT(2),XT(3) ) )
YTMIN = MIN(YT(1),MIN( YT(2),YT(3) ) )
IF (XP .GE. XTMIN .AND. XP .LE. XTMAX .AND. &
YP .GE. YTMIN .AND. YP .LE. YTMAX) THEN
CALL PINPOK(XP,YP,3,XT,YT,INTRI)
IF (INTRI .EQ. 1) THEN
NRFIND = K
NROLDFIND = NRFIND
do idim=1,NDIM
ZT(idim,1) = ZS(idim,INDX(1,K))
ZT(idim,2) = ZS(idim,INDX(2,K))
ZT(idim,3) = ZS(idim,INDX(3,K))
end do
CALL LINEAR (XT, YT, ZT, NDIM, XP, YP, ZP, JSLO, SLO, JATEK, wf)
ind(1) = ik1
ind(2) = ik2
ind(3) = ik3
RETURN
ENDIF
ENDIF
ENDDO
if (k1 == 1 .and. k2 == numtri) then
NROLDFIND = numtri / 2
return
endif
IF (NRfind == 0) THEN
GOTO 5
ENDIF
RETURN
END
!> find triangle for interpolation with kdtree
!> will initialize kdtree and triangulation connectivity
subroutine findtri_kdtree(XP,YP,ZP,XS,YS,ZS,NS,NDIM,NRFIND,INTRI,JSLO,SLO,JATEK,jadum,ierror,ind,wf)
use m_missing
use m_triangle
use m_kdtree2
use MessageHandling
use m_sferic
use m_alloc
implicit none
double precision, intent(in) :: xp, yp !< node coordinates
double precision, intent(out) :: zp(NDIM) !< node values
integer, intent(in) :: NS !< number of samples
integer, intent(in) :: NDIM !< sample vector dimension
double precision, intent(in) :: xs(NS), ys(NS), zs(NDIM,NS) !< sample coordinates and values
integer, intent(out) :: NRFIND !< triangle index
integer, intent(out) :: intri !< in triangle (1) or not (0)
integer, intent(in) :: jslo !< get slope (1) or not (0)
double precision, intent(out) :: slo(NDIM) !< slope
integer, intent(in) :: JATEK !< draw to screen (1) or not (0)
integer, intent(in) :: jadum !< debug (1) or not (0)
integer, intent(out) :: ierror !< error (1) or not (0)
integer, intent(out) :: ind(3) !<
double precision, intent(out) :: wf(3) !<
double precision, dimension(:), allocatable :: xx, yy !< triangle circumcenter coordinates, dim(numtri)
double precision, dimension(3) :: xv, yv ! triangle node coordinates
double precision, dimension(NDIM,3) :: zv ! triangle node vectors
double precision :: xz, yz ! triangle circumcenter
double precision :: SL,SM,XCR,YCR,CRP
integer :: NN ! number of nearest points
integer :: IDIM ! dimensionality
integer :: i, inod, inod1, inod2, inod3, inod4, ii, iii, iip1, inext, j, k
integer :: iothertri, iiothertri, inodother, numtrisnew, jacros, iothertriangle
integer :: jsferic_store, ierr
integer :: iedge, k1, k2, numsearched
double precision, parameter :: dtol = 1d-8
integer, parameter :: MAXTRICON=100
integer, parameter :: INIT_TRIS_PER_NODE=6
double precision, external :: dcosphi
ierror = 1
NRFIND = 0
if ( numtri.le.1 ) then
! call qnerror('findtri_kdtree: numtri<=1', ' ', ' ')
goto 1234
end if
! store jsferic
jsferic_store = jsferic
if ( treeglob%itreestat /= ITREE_OK ) then
! INITIALIZE kdetree2
call mess(LEVEL_INFO,'Determining triangle connectivity and computing polygons and circumcenters...')
allocate(xx(numtri), yy(numtri), stat=ierr)
call aerr('xx(numtri), yy(numtri)', ierr, 2*numtri)
allocate(imask(numtri), stat=ierr)
call aerr('imask(numtri)', ierr, numtri)
imask = 0
IDENT = 0
allocate(LNtri(2,numedge), stat=ierr)
call aerr('LNtri(2,numedge)', ierr, 2*numedge)
LNtri = 0
! dlaun is not spherical proof: set global jsferic for circumcenter
jsferic = 0
do i=1,numtri
do ii=1,3
! generate edge-triangle connectivity
iedge = triedge(ii,i)
if ( LNtri(1,iedge).eq.0 ) then
LNtri(1,iedge) = i
else
LNtri(2,iedge) = i
end if
inod = indx(ii,i)
! get triangle polygon
xv(ii) = xs(inod)
yv(ii) = ys(inod)
zv(1:NDIM,ii) = zs(1:NDIM,inod)
end do
! compute triangle circumcenter
!call circumcenter3(3, xv, yv, xx(i), yy(i))
xx(i) = sum(xv(1:3)) / 3d0
yy(i) = sum(yv(1:3)) / 3d0
end do
! restore jsferic
jsferic = jsferic_store
call mess(LEVEL_INFO, 'done')
call build_kdtree(treeglob, numtri, xx, yy, ierror)
! deallocate
deallocate(xx, yy)
if ( ierror.ne.0 ) then
goto 1234
end if
end if
if ( jadum.eq.1 ) then
continue
end if
! update identifier
IDENT=IDENT+1
! fill query vector
call make_queryvector_kdtree(treeglob,xp,yp)
! get first triangle
call kdtree2_n_nearest(treeglob%tree,treeglob%qv,1,treeglob%results)
inext = treeglob%results(1)%idx
! perform a line search
numsearched = 1
i = 0
mainloop:do while ( NRFIND.eq.0 .and. numsearched.le.2*NUMTRI .and. numsearched.gt.0 )
inext = treeglob%results(1)%idx
numsearched = 0
! check current triangle
do ii=1,3
inod = indx(ii,inext)
xv(ii) = xs(inod)
yv(ii) = ys(inod)
zv(1:NDIM,ii) = zs(1:NDIM,inod)
end do
! get cell centroid
xz = sum(xv(1:3)) / 3d0
yz = sum(yv(1:3)) / 3d0
do while ( inext.gt.0 .and. inext.le.numtri .and. numsearched.le.2*NUMTRI ) ! numsearched: safety
i = inext
! check current triangle
do ii=1,3
inod = indx(ii,i)
xv(ii) = xs(inod)
yv(ii) = ys(inod)
zv(1:NDIM,ii) = zs(1:NDIM,inod)
end do
if ( jadum.eq.1 ) then
call tektri(xv,yv,31)
xz = sum(xv(1:3)) / 3d0
yz = sum(yv(1:3)) / 3d0
call cirr(xz,yz,31)
end if
if ( imask(i).eq.IDENT ) then
intri = 0
else
numsearched = numsearched+1
call pinpok(xp,yp,3,xv,yv,intri)
imask(i) = IDENT
end if
if ( intri.eq.1 ) then
NRFIND = i
exit mainloop
end if
! dlaun is not spherical proof: set global jsferic
jsferic = 0
! proceed to next triangle, which is adjacent to the edge that is cut by the line from the current triangle to the query point
xz = sum(xv(1:3)) / 3d0
yz = sum(yv(1:3)) / 3d0
inext = 0
if ( jadum.eq.1 ) then
call setcol(221)
call movabs(xz,yz)
call lnabs(xp,yp)
call setcol(31)
end if
do ii=1,3
iedge=triedge(ii,i)
if ( LNtri(2,iedge).eq.0 ) cycle
iothertriangle = LNtri(1,iedge) + LNtri(2,iedge) - i
if ( imask(iothertriangle).eq.IDENT ) then
cycle
end if
k1 = edgeindx(1,iedge)
k2 = edgeindx(2,iedge)
call CROSS(xz, yz, xp, yp, xs(k1), ys(k1), xs(k2), ys(k2), JACROS,SL,SM,XCR,YCR,CRP)
! use tolerance
if ( jacros.eq.0 ) then
IF (SM .GE. 0d0-dtol .AND. SM .LE. 1d0+dtol .AND. &
SL .GE. 0d0-dtol .AND. SL .LE. 1d0+dtol) THEN
JACROS = 1
ENDIF
end if
if ( jadum.eq.1 .and. jacros.eq.1 ) then
call movabs(xs(k1),ys(k1))
call lnabs(xs(k2),ys(k2))
call qnerror(' ', ' ', ' ')
end if
if ( jacros.eq.1 ) then ! only select this edge if it has a second adjacent triangle
inext = iothertriangle
exit
end if
end do
! restore jsferic
jsferic = jsferic_store
!if ( jadum.eq.1 ) then
! call tektri(xv,yv,0)
!end if
end do
end do mainloop
if ( intri.eq.0 .and. inext.gt.0 ) then
if ( numsearched.gt.10 ) write(6,"('numsearched= ', I0)") numsearched
! call qnerror('findtri_kdtree: error', ' ', ' ')
end if
if (intri .eq. 1) then
call linear(xv, yv, zv, NDIM, xp, yp, zp, JSLO, SLO, JATEK, wf)
do k = 1,3
ind(k) = indx(k,nrfind)
enddo
endif
ierror = 0
1234 continue
return
end subroutine findtri_kdtree
SUBROUTINE LINEAR ( X, Y, Z, NDIM, XP, YP, ZP, JSLO, SLO, JATEK, wf)
USE M_MISSING
! use unstruc_colors
implicit none
double precision :: a11
double precision :: a12
double precision :: a21
double precision :: a22
double precision :: a31
double precision :: a32
double precision :: b1
double precision :: b2
double precision :: det
double precision :: dum
integer :: jatek
integer :: jslo
integer :: ncol
integer :: ndraw
double precision :: r3
double precision :: rlam
double precision :: rmhu
double precision :: x1
double precision :: x2
double precision :: x3
double precision :: xmax
double precision :: xmin
double precision :: xn
double precision :: xp
double precision :: xy
double precision :: y1
double precision :: y2
double precision :: y3
double precision :: ymax
double precision :: ymin
double precision :: yn
double precision :: yp
double precision :: z3
double precision :: zn
integer :: idim
integer :: NDIM !< sample vector dimension
double precision :: X(3),Y(3),Z(NDIM,3), wf(3)
double precision :: zp(NDIM)
double precision :: slo(NDIM)
COMMON /DRAWTHIS/ NDRAW(40)
double precision :: getdx, getdy
ZP = dmiss
A11 = getdx(x(1),y(1),x(2),y(2)) ! X(2) - X(1)
A21 = getdy(x(1),y(1),x(2),y(2)) ! Y(2) - Y(1)
A12 = getdx(x(1),y(1),x(3),y(3)) ! X(3) - X(1)
A22 = getdy(x(1),y(1),x(3),y(3)) ! Y(3) - Y(1)
B1 = getdx(x(1),y(1),xp ,yp ) ! XP - X(1)
B2 = getdy(x(1),y(1),xp ,yp ) ! YP - Y(1)
DET = A11 * A22 - A12 * A21
IF (ABS(DET) .LT. 1E-12) THEN ! Jan Mooiman 07-01-2015
RETURN
ENDIF
RLAM = ( A22 * B1 - A12 * B2) / DET
RMHU = (-A21 * B1 + A11 * B2) / DET
wf(3) = rmhu
wf(2) = rlam
wf(1) = 1d0 - rlam - rmhu
ZP = Z(:,1) + RLAM * (Z(:,2) - Z(:,1)) + RMHU * (Z(:,3) - Z(:,1))
IF (JATEK .EQ. 1) THEN
CALL ISOCOL(ZP(1),NCOL)
CALL KCIR(XP,YP,ZP(1))
CALL TEKTRI(X,Y,NCOL)
IF (MAX(ABS(A21),ABS(A22)) .GT. 500) THEN
DUM = 0
ENDIF
ENDIF
IF (JSLO .EQ. 1) THEN
do idim = 1,NDIM
A31 = Z(idim,2) - Z(idim,1)
A32 = Z(idim,3) - Z(idim,1)
X3 = (A21*A32 - A22*A31)
Y3 = -(A11*A32 - A12*A31)
Z3 = (A11*A22 - A12*A21)
R3 = SQRT(X3*X3 + Y3*Y3 + Z3*Z3)
IF (R3 .NE. 0) THEN
XN = X3/R3
YN = Y3/R3
ZN = Z3/R3
XY = SQRT(XN*XN + YN*YN)
IF (ZN .NE. 0) THEN
SLO(idim) = ABS(XY/ZN)
ELSE
SLO(idim) = dmiss
ENDIF
ELSE
SLO(idim) = dmiss
ENDIF
end do
ENDIF
RETURN
END
subroutine read_samples_from_dem(filnam, jadoorladen)
use dem
use m_missing
use m_samples
implicit none
character(len=*), intent(in) :: filnam
integer, intent(in) :: jadoorladen
integer :: ndraw
COMMON /DRAWTHIS/ NDRAW(40)
integer :: i, j, is, istep
type(DEMInfo) :: dem_info
integer, allocatable :: arr(:,:)
double precision, allocatable :: xarr(:,:), yarr(:,:)
character(len=10) :: TEX
call savesam()
if (jadoorladen == 0) then
ns = 0
end if
call read_dem_file(trim(filnam),dem_info, xarr, yarr, arr)
if (dem_info%rows <= 0 .or. dem_info%cols <= 0) then
call message('No samples read from file ', filnam, ' ')
return
end if
call increasesam(ns + dem_info%rows*dem_info%cols)
WRITE(TEX,'(I10)') dem_info%rows*dem_info%cols
CALL READYY('Filtering '//TRIM(TEX)//' Samples Points',0d0)
istep = int(dem_info%rows/100d0)
do i=1,dem_info%rows
do j=1,dem_info%cols
if (arr(i,j) == NODATA) then
continue
else
ns = ns + 1
xs(ns) = xarr(i,j)
ys(ns) = yarr(i,j)
zs(ns) = dble(arr(i,j))
end if
end do
IF (MOD(i,istep) .EQ. 0) THEN
CALL READYY(' ',MIN( 1d0,dble(i)/dem_info%rows) )
ENDIF
end do
deallocate(xarr, yarr, arr)
CALL READYY(' ',-1d0)
IF (NS .GT. 100000) NDRAW(32) = 7 ! Squares (faster than circles)
IF (NS .GT. 500000) NDRAW(32) = 3 ! Small dots (fastest)
WRITE(TEX,'(I10)') NS
CALL READYY('Sorting '//TRIM(TEX)//' Samples Points',0d0)
IF (NS .GT. 1) THEN
CALL TIDYSAMPLES(XS,YS,ZS,IPSAM,NS,MXSAM,MYSAM)
call get_samples_boundingbox()
IPSTAT = IPSTAT_OK
END IF
CALL READYY(' ',-1d0)
end subroutine read_samples_from_dem
subroutine read_samples_from_arcinfo(filnam, jadoorladen)
use m_missing
use m_samples
use m_samples_refine, only: iHesstat, iHesstat_DIRTY
use m_arcinfo
use unstruc_display, only: jagui
implicit none
character(len=*), intent(in) :: filnam
integer, intent(in) :: jadoorladen
integer :: i, j, istep, marc
character(len=10) :: TEX
integer :: ndraw
COMMON /DRAWTHIS/ NDRAW(40)
CALL READYY('Reading arcinfo file',0d0)
call oldfil(marc, filnam)
call reaarc (marc,jagui)
CALL DOCLOSE(marc)
CALL READYY('Reading arcinfo file',1d0)
if (mca <= 0 .or. nca <= 0) then
call message('No samples read from file ', filnam, ' ')
return
end if
call savesam()
if (jadoorladen == 0) then
ns = 0
end if
CALL INCREASEsam(ns+mca*nca)
WRITE(TEX,'(I10)') mca*nca
CALL READYY('Filtering '//TRIM(TEX)//' Samples Points',0d0)
istep = max(int(mca/100d0+.5d0),1)
! SPvdP: j needs to be fastest running index
do i = 1,mca
IF (MOD(i,istep) .EQ. 0) THEN
CALL READYY('Filtering '//TRIM(TEX)//' Samples Points',min( 1d0,dble(i)/mca))
ENDIF
do j = nca,1,-1 ! SPvdP: first line needs to be nca'th row
! if (d(I,J) .ne. dmiss) then ! SPvdP: we need to maintain structured data
ns = ns+1
xs(ns) = x0 + dxa*(i-1)
ys(ns) = y0 + dya*(j-1)
zs(ns) = d(i,j)
! endif
enddo
enddo
CALL READYY(' ',-1d0)
! mark samples as structured, and in supply block sizes
MXSAM = nca ! j is fastest running index
MYSAM = mca
IPSTAT = IPSTAT_NOTOK
! new sample set: no Hessians computed yet
iHesstat = iHesstat_DIRTY
deallocate(d) ! Save memory, arcinfo block is no longer needed.
IF (NS .GT. 100000) NDRAW(32) = 7 ! Squares (faster than circles)
IF (NS .GT. 500000) NDRAW(32) = 3 ! Small dots (fastest)
! No TIDYSAMPLES required: arcinfo grid was already loaded in correctly sorted order.
do i=1,NS
IPSAM(i) = i
end do
call get_samples_boundingbox()
IPSTAT = IPSTAT_OK
end subroutine read_samples_from_arcinfo
SUBROUTINE REASAM(MSAM, JADOORLADEN)
USE M_MISSING
USE M_SAMPLES
use m_alloc
implicit none
integer, intent(inout) :: msam !< already opened file pointer to sample file
integer, intent(in) :: jadoorladen !< whether to append to global set (1) or start empty (0)
integer :: ierr
integer :: jflow
integer :: jqn
integer :: mcs
integer :: ncs
integer :: ndraw
integer :: nkol
integer :: nrow
integer :: ns1
integer :: nsm
integer :: num
integer :: K, K0
double precision :: x, y, z
double precision :: XX, YY, ZZ, ZZ2
COMMON /PHAROSFLOW/ JFLOW
COMMON /PHAROSLINE/ REC1
COMMON /SAMPLESADM/ MCS,NCS,NS1
COMMON /QNRGF/ JQN
COMMON /DRAWTHIS/ NDRAW(40)
CHARACTER REC*132, TEX*10, REC1*132
LOGICAL THISISANUMBER
CALL SAVESAM()
NSM = 0
MXSAM = 0
MYSAM = 0
IPSTAT = IPSTAT_NOTOK
nkol = 0
CALL READYY('Counting nr. of Samples ',0d0)
11 READ (MSAM,'()',END = 31)
NSM = NSM + 1
GOTO 11
31 NSMAX = 1.2d0*(NSM + JADOORLADEN*NS)
IF (NSMAX .GT. 100000) NDRAW(32) = 7
IF (NSMAX .GT. 500000) NDRAW(32) = 3
IF (ALLOCATED (XS) ) DEALLOCATE (XS,YS,ZS)
ALLOCATE (XS(NSMAX),YS(NSMAX),ZS(NSMAX),STAT=IERR)
CALL AERR ('XS(NSMAX),YS(NSMAX),ZS(NSMAX)',IERR,NSMAX)
if ( allocated(ipsam) ) deallocate(ipsam)
allocate(ipsam(NSMAX),stat=ierr)
call aerr('ipsam(NSMAX)',ierr,NSMAX)
CALL READYY(' ',-1d0)
REWIND(MSAM)
WRITE(TEX,'(I10)') NSM
CALL READYY('Reading '//TRIM(TEX)//' Sample Points',0d0)
IF (JADOORLADEN .EQ. 0) THEN
CALL XMISAR(XS,NSMAX)
CALL XMISAR(YS,NSMAX)
CALL MISAR(ZS,NSMAX)
K = 0
ELSE
CALL RESTORESAM()
K = NS
NS1 = NS
ENDIF
K0 = K
! check of dit een PHAROS file is
JFLOW = 1
14 READ (MSAM,'(A)',END = 30) REC1
if (rec1(1:1) == '*') goto 14
IF ( .NOT. (THISISANUMBER(REC1)) ) THEN
READ (MSAM,'(A)',END = 30) REC
IF ( THISISANUMBER(REC) ) THEN
READ (REC,*,ERR = 16) NROW,NKOL
GOTO 15
16 CONTINUE
READ (MSAM,'(A)',END = 30) REC
READ (REC,*,ERR = 15) NUM,X,Y,Z
JFLOW = 3
ENDIF
ENDIF
15 CONTINUE
REWIND (MSAM)
KMOD = max(1, NSM/100)
10 CONTINUE
READ (MSAM,'(A)',END = 30) REC
IF (REC(1:1) .EQ. '*') GOTO 10
IF ( .NOT. (THISISANUMBER(REC)) ) THEN
! we nemen aan dat er net een blokcode is gelezen
! en we lezen meteen de nrow ncol regel, maar checken die regel niet
READ (MSAM,'(A)',END = 30) REC
ELSE
IF (JFLOW .EQ. 3) THEN
READ (REC,*,ERR = 40) NUM,XX,YY,ZZ
ELSE IF (NKOL == 4) THEN
READ (REC,*,ERR = 40) XX,YY, ZZ, ZZ2
if (zz .ne. -999d0) then
zz = sqrt(zz*zz + zz2*zz2)
endif
ELSE
READ (REC,*,end = 40) XX,YY,ZZ
READ (REC,*,ERR = 40) XX,YY,ZZ
ENDIF
IF (K .LE. NSMAX-1 .AND. XX .NE. XYMIS .AND. &
ZZ .NE. dmiss .AND. ZZ .NE. 999.999d0 .and. &
.not.(isnan(XX) .or. isnan(YY) .or. isnan(ZZ)) ) THEN
K = K + 1
NS = K
XS(K) = XX
YS(K) = YY
ZS(K) = ZZ
ENDIF
IF (MOD(K-K0,KMOD) .EQ. 0) THEN
CALL READYY(' ',MIN( 1d0,dble(K)/NSM) )
ENDIF
ENDIF
GOTO 10
40 CONTINUE
WRITE(TEX,'(I10)') K
CALL QNERROR('ERROR READING SAMPLES FILE LINE NR ',TEX,REC)
30 CONTINUE
IF (K .GT. NSMAX) THEN
WRITE(TEX,'(I8)') NSMAX
CALL QNERROR('ONLY',TEX,'SAMPLE POINTS CAN BE LOADED')
WRITE(TEX,'(I8)') K
CALL QNERROR('YOU TRIED TO LOAD',TEX,'SAMPLE POINTS')
ENDIF
CALL READYY(' ',-1d0)
WRITE(TEX,'(I10)') NS
CALL READYY('Sorting '//TRIM(TEX)//' Samples Points',0d0)
IF (NS .GT. 1) THEN
CALL TIDYSAMPLES(XS,YS,ZS,IPSAM,NS,MXSAM,MYSAM)
call get_samples_boundingbox()
IPSTAT = IPSTAT_OK
END IF
CALL READYY(' ',-1d0)
call doclose(MSAM)
RETURN
END
SUBROUTINE INCREASESAM(N)
USE M_SAMPLES
USE M_MISSING
use m_alloc
implicit none
integer, intent(in) :: n !< New size for sample set #3.
integer :: ierr
IF (N < NSMAX) RETURN
NSMAX = MAX(10000,INT(1.2d0*N))
call realloc(xs, NSMAX, keepExisting=.true., fill = dmiss, stat=ierr)
call realloc(ys, NSMAX, keepExisting=.true., fill = dmiss, stat=ierr)
call realloc(zs, NSMAX, keepExisting=.true., fill = dmiss, stat=ierr)
call realloc(ipsam, NSMAX, keepExisting=.false., fill=0, stat=ierr)
CALL AERR ('XS(NSMAX),YS(NSMAX),ZS(NSMAX),IPSAM(NSMAX)',IERR,NSMAX)
! user is editing samples: mark samples as unstructured
MXSAM = 0
MYSAM = 0
IPSTAT = IPSTAT_NOTOK
END
SUBROUTINE INCREASESAM3(N)
USE M_SAMPLES3
USE M_MISSING
use m_alloc
implicit none
integer, intent(in) :: n !< New size for sample set #3.
integer :: ierr
integer :: nsmax
NSMAX = SIZE(XS3)
IF (N < NSMAX) RETURN
NSMAX = MAX(10000,INT(1.2d0*N))
call realloc(xs3, NSMAX, keepExisting=.true., fill = dmiss, stat=ierr)
call realloc(ys3, NSMAX, keepExisting=.true., fill = dmiss, stat=ierr)
call realloc(zs3, NSMAX, keepExisting=.true., fill = dmiss, stat=ierr)
CALL AERR ('XS3(NSMAX),YS3(NSMAX),ZS3(NSMAX)',IERR,NSMAX)
END
SUBROUTINE SAVESAM()
USE M_SAMPLES2
USE M_SAMPLES
USE M_MISSING
use m_alloc
implicit none
integer :: ierr
NS2 = NS
MXSAM2 = MXSAM
MYSAM2 = MYSAM
IF (NS .EQ. 0) RETURN
call realloc(xs2, ns, keepExisting=.false., fill=dmiss, stat=ierr)
call realloc(ys2, ns, keepExisting=.false., fill=dmiss, stat=ierr)
call realloc(zs2, ns, keepExisting=.false., fill=dmiss, stat=ierr)
CALL PUTAR(XS,XS2,NS)
CALL PUTAR(YS,YS2,NS)
CALL PUTAR(ZS,ZS2,NS)
! NS2 = NS
RETURN
END
SUBROUTINE RESTORESAM()
USE M_SAMPLES2
USE M_SAMPLES
implicit none
MXSAM = 0 ! unstructured samples by default
MYSAM = 0
IPSTAT = IPSTAT_NOTOK
NS = NS2
MXSAM = MXSAM2
MYSAM = MYSAM2
IF (NS2 == 0) RETURN
CALL PUTAR(XS2,XS,NS2)
CALL PUTAR(YS2,YS,NS2)
CALL PUTAR(ZS2,ZS,NS2)
RETURN
END
!> Increase size of global polyline array.
!! Specify new size and whether existing points need to be maintained.
SUBROUTINE INCREASEPOL(N, jaKeepExisting)
USE M_POLYGON
USE M_MISSING
use m_alloc
implicit none
integer :: n !< Desired new minimum size
integer :: jaKeepExisting !< Whether or not (1/0) to keep existing points.
logical :: jakeep
integer :: maxpolcur
integer :: ierr
maxpolcur = size(xpl)
IF (N < maxpolcur ) THEN
RETURN
ENDIF
MAXPOL = MAX(100000,INT(5d0*N))
jakeep = jaKeepExisting==1
call realloc(xpl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
call realloc(ypl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
call realloc(zpl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
if (jakol45 == 1) then
call realloc(dzl, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
call realloc(dzr, maxpol, keepExisting=jakeep, fill=dxymis, stat=ierr)
endif
! make sure nampli is allocated
if ( .not.allocated(nampli) ) then
allocate(nampli(0))
end if
end subroutine increasepol
!> Copies the global polygon into the backup polygon arrays.
SUBROUTINE SAVEPOL()
USE M_POLYGON
use m_alloc
use m_missing
implicit none
call realloc(xph, maxpol, keepExisting=.false.)
call realloc(yph, maxpol, keepExisting=.false.)
call realloc(zph, maxpol, keepExisting=.false.)
IF (NPL > 0) THEN
XPH(1:NPL) = XPL(1:NPL)
YPH(1:NPL) = YPL(1:NPL)
ZPH(1:NPL) = ZPL(1:NPL)
ENDIF
MPS = MP
NPH = NPL
RETURN
END subroutine savepol
!> Puts back a previously saved backup polygon into the global polygon arrays.
SUBROUTINE RESTOREPOL()
USE M_POLYGON
use m_alloc
use m_missing
implicit none
maxpol = max(maxpol, nph)
call realloc(xpl, maxpol, keepExisting=.false.)
call realloc(ypl, maxpol, keepExisting=.false.)
call realloc(zpl, maxpol, keepExisting=.false.)
IF (NPH > 0) THEN
XPL(1:NPH) = XPH(1:NPH)
YPL(1:NPH) = YPH(1:NPH)
ZPL(1:NPH) = ZPH(1:NPH)
ENDIF
MP = MPS
NPL = NPH
RETURN
end subroutine restorepol
SUBROUTINE XMISAR (X, MMAX)
USE M_MISSING
implicit none
integer :: i
integer :: mmax
double precision :: x
DIMENSION X(MMAX)
DO 10 I = 1,MMAX
X(I) = XYMIS
10 CONTINUE
RETURN
END
SUBROUTINE TIDYSAMPLES(XS,YS,ZS,IPSAM,NS,MXSAM,MYSAM)
implicit none
integer :: ns
double precision :: XS(NS), YS(NS), ZS(NS) !< sample coordinates
integer, dimension(NS), intent(out) :: IPSAM !< permutation array (increasing x-coordinate)
integer, intent(in) :: MXSAM, MYSAM !< structured sample data dimensions (>0) or unstructured (0)
! IF (NS .GT. 1) CALL RSORT3(XS,YS,ZS,NS)
if ( NS.gt.1 ) then
call indexx(Ns,xs,IPSAM)
end if
! remove double/missing samples (non-structured sample data only)
if ( MXSAM*MYSAM.ne.NS ) then
CALL READYY(' ',0.3d0)
IF (NS .GT. 1) CALL RMDOUBLE(XS,YS,ZS,IPSAM,NS)
end if
CALL READYY(' ',1d0)
RETURN
END
!> determine sample bounding box
subroutine get_samples_boundingbox()
use m_samples
use m_missing
implicit none
integer :: i
xsammin = huge(1d0)
xsammax = -huge(1d0)
ysammin = huge(1d0)
ysammax = -huge(1d0)
do i=1,NS
if ( xs(i).ne.DMISS .and. ys(i).ne.DMISS .and. zs(i).ne.DMISS ) then
xsammin = min(xsammin,xs(i))
xsammax = max(xsammax,xs(i))
ysammin = min(ysammin,ys(i))
ysammax = max(ysammax,ys(i))
end if
end do
return
end subroutine get_samples_boundingbox
SUBROUTINE RSORT3new (X, Y, Z, N) ! 1 !! second faster than
implicit none
double precision :: X(N), Y(N), Z(N)
integer, allocatable :: ind(:)
double precision, allocatable :: h(:)
integer :: k, n
allocate(ind(n), h(n))
call indexx(n,x,ind)
h = x
do k = 1,n
x(k) = h(ind(k))
enddo
h = y
do k = 1,n
y(k) = h(ind(k))
enddo
h = z
do k = 1,n
z(k) = h(ind(k))
enddo
deallocate(ind,h)
end SUBROUTINE RSORT3new
SUBROUTINE RSORT3(X, Y, Z, N)
implicit none
integer :: j
integer :: j1
integer :: k0
integer :: kk
integer :: lk
integer :: n
integer :: nm
double precision :: temp
double precision :: X(N), Y(N), Z(N)
IF (N .EQ. 0) RETURN
LK = N / 2
K0 = LK
KK = K0
20 J = 2 * KK
J1 = J + 1
30 IF ( J1 .LE. N ) THEN
IF ( X(J) .LT. X(J1) ) J = J1
ENDIF
IF ( X(KK) .LT. X(J) ) THEN
TEMP = X(J)
X(J) = X(KK)
X(KK) = TEMP
TEMP = Y(J)
Y(J) = Y(KK)
Y(KK) = TEMP
TEMP = Z(J)
Z(J) = Z(KK)
Z(KK) = TEMP
IF ( J .LE. LK ) THEN
KK = J
GOTO 20
ENDIF
ENDIF
K0 = K0 - 1
IF ( K0 .NE. 0 ) THEN
KK = K0
J = 2 * KK
J1 = J + 1
GOTO 30
ENDIF
NM = N
65 TEMP = X(1)
X(1) = X(NM)
X(NM) = TEMP
TEMP = Y(1)
Y(1) = Y(NM)
Y(NM) = TEMP
TEMP = Z(1)
Z(1) = Z(NM)
Z(NM) = TEMP
NM = NM - 1
IF ( NM .EQ. 1 ) RETURN
KK = 1
70 J = 2 * KK
J1 = J + 1
IF ( J .GT. NM ) GOTO 65
IF ( J1 .LE. NM .AND. X(J) .LT. X(J1) ) J = J1
IF ( X(KK) .GE. X(J) ) GOTO 65
TEMP = X(J)
X(J) = X(KK)
X(KK) = TEMP
TEMP = Y(J)
Y(J) = Y(KK)
Y(KK) = TEMP
TEMP = Z(J)
Z(J) = Z(KK)
Z(KK) = TEMP
KK = J
GOTO 70
END
SUBROUTINE RMDOUBLE(XS,YS,ZS,IPSAM,NS)
USE M_MISSING
use m_kdtree2
use m_sferic
use unstruc_messages
implicit none
integer :: i
integer :: j
integer :: jadouble
integer :: k
integer :: ns
integer :: nsorg
integer :: numweg
integer :: isam, jsam, ksam
double precision, intent(inout) :: XS(NS), YS(NS), ZS(NS)
integer, dimension(NS), intent(inout) :: IPSAM !< permutation array (increasing x-coordinate)
integer, dimension(:), allocatable :: newnode
double precision, dimension(:), allocatable :: xx, yy ! non-missing sample coordinates
integer, dimension(:), allocatable :: iperm ! permutation array
double precision :: t0, t1, t2, t3, t4
integer :: ii, jj, NN, num, nummerged, ierror
integer :: jakdtree=1
integer :: jsferic_store
character(len=128) :: txt
double precision, parameter :: dtol2 = 1d-8 ! sample-on-top of each other tolerance, squared
CHARACTER OUD*8, NIEUW*8
NSORG = NS
allocate(newnode(NS))
! store jsferic
jsferic_store = jsferic
! safety
!if ( NS.lt.100 ) then
! jakdtree=0
!end if
call klok(t0)
if ( jakdtree.eq.1 ) then
! force Cartesian coordinates
jsferic = 0
! get non-missing sample coordinates
allocate(xx(NS))
xx = 0d0
allocate(yy(NS))
yy = 0d0
allocate(iperm(NS))
iperm = 0
! get non-missing sample coordinates
num = 0
do i=1,NS
if ( xs(i).ne.DMISS .and. ys(i).ne.DMISS ) then
! fix for Wim: already check samples with latest (not understood kdtree error)
if ( num.gt.1 ) then
if( xs(i).eq.xs(num) .and. ys(i).eq.ys(num) ) cycle
end if
num = num+1
xx(num) = xs(i)
yy(num) = ys(i)
iperm(num) = i
end if
end do
! initialize kdtree
call build_kdtree(treeglob,num,xx,yy,ierror)
! deallocate arrays with non-missing node coordinates
deallocate(xx)
deallocate(yy)
if ( ierror.ne.0 ) then
! deallocate permutation array
if ( allocated(iperm) ) deallocate(iperm)
! deallocate kdtree
if ( treeglob%itreestat.ne.ITREE_EMPTY ) call delete_kdtree2(treeglob)
! disable kdtree
jakdtree = 0
end if
end if
nummerged=0
! fill double samples with DMISS
5 CONTINUE
JADOUBLE = 0
if ( jakdtree.eq.1 ) then
! find samples on top of each other
do ii=1,num
i=iperm(ii)
if ( i.eq.0 ) cycle ! already merged
! fill query vector
call make_queryvector_kdtree(treeglob,xs(i),ys(i))
! count number of points in search area
NN = kdtree2_r_count(treeglob%tree,treeglob%qv,dtol2)
if ( NN.gt.1 ) then ! at least two samples need to be merged
! resize results array if necessary
call realloc_results_kdtree(treeglob,NN)
! find other nodes
call kdtree2_n_nearest(treeglob%tree,treeglob%qv,NN,treeglob%results)
! merge with other nodes
do k=1,NN
jj = treeglob%results(k)%idx
j = iperm(jj)
! exclude own sample and samples already deleted
if ( j.ne.i .and. j.gt.0 ) then
if ( xs(i).eq.xs(j) .and. ys(i).eq.ys(j) ) then ! NOT SPHERICAL-PROOF
iperm(jj) = 0
xs(j) = DMISS
jadouble = 1
nummerged = nummerged+1
end if
end if
end do
end if
end do
else ! non kdtree
DO 10 I = 1,NS-1
ISAM = IPSAM(I)
IF (XS(ISAM) .NE. XYMIS .and. ZS(ISAM) .NE. DMISS) THEN
J = I
15 CONTINUE
IF (J .LT. NS) THEN
J = J + 1
JSAM = IPSAM(J)
IF (XS(ISAM) .EQ. XS(JSAM)) THEN
IF (YS(ISAM) .EQ. YS(JSAM)) THEN
XS(JSAM) = XYMIS
JADOUBLE = 1
ENDIF
GOTO 15
ENDIF
ENDIF
ENDIF
10 CONTINUE
end if
call klok(t1)
! remove double samples
K = 0
newnode = 0
DO 20 I = 1,NS
IF (XS(I) .NE. XYMIS .and. ZS(I) .NE. DMISS) THEN
K = K + 1
XS(K) = XS(I)
YS(K) = YS(I)
ZS(K) = ZS(I)
newnode(i) = k ! old-to-new sample number
ENDIF
20 CONTINUE
call klok(t2)
! update permutation array
k = 0
do i=1,NS
j = IPSAM(i) ! old node number
if ( newnode(j).gt.0 ) then
k=k+1
IPSAM(k) = newnode(j)
end if
end do
call klok(t3)
! set new number of samples
NS = K
IF (JADOUBLE .EQ. 1) GOTO 5
NUMWEG = NSORG - K
IF (NUMWEG .GE. 1) THEN
WRITE(OUD,'(I8)') NUMWEG
! CALL QNERROR('NUMBER OF DOUBLE POINTS REMOVED',OUD,' ')
ENDIF
call klok(t4)
! output message
if ( jakdtree.eq.1 ) then
txt = ''
write(txt, "('merged ', I0, ' samples in ', F0.2, ' seconds.')") nummerged, t4-t0
call mess(LEVEL_INFO, trim(txt))
end if
1234 continue
! deallocate
if ( allocated(newnode) ) deallocate(newnode)
if ( jakdtree.eq.1 ) then
! deallocate permutation array
if ( allocated(iperm) ) deallocate(iperm)
! deallocate kdtree
if ( treeglob%itreestat.ne.ITREE_EMPTY ) call delete_kdtree2(treeglob)
end if
! restore jsferic
jsferic = jsferic_store
RETURN
END
SUBROUTINE DLAUN(XS,YS,NS,jatri,ierr)
USE M_TRIANGLE
use m_sferic
use m_alloc
use unstruc_messages
implicit none
integer, intent(out) :: ierr
integer :: maxtri
integer :: nh
integer :: ns
integer, intent(in) :: jatri !< Type of DLaun triangulation: 1: just triangulate,
!! 3: also produce node-edge-triangle mapping tables
!! for use in Triangulatesamplestonetwork.
integer :: nsm
double precision :: trisize
PARAMETER (NH = 1) ! SPvdP: too small if jatri.eq.0
double precision :: XS(ns), YS(ns)
double precision :: XH(NH), YH(NH)
integer, allocatable, dimension(:) :: idum
call mess(LEVEL_INFO,'Starting Delaunay triangulation...')
! check memory
allocate ( idum(50*Ns) ,stat=ierr) ! probably not enough
call aerr('idum(50*Ns)',ierr,-50*Ns)
if ( ierr.ne.0 ) then
call qnerror('dlaun: out of memory', ' ', ' ')
! if ( allocated(idum) ) deallocate(idum) ! gives an error
return
end if
deallocate(idum)
NUMTRI = 0
NUMEDGE = 0
IF (NS .LT. 3) RETURN
if (jatri /= 1 .and. jatri /= 3) then
call mess(LEVEL_INFO, 'Invalid jatri value in DLAUN:', jatri)
return
end if
numtri = -1
do while ( numtri.lt.0 )
NSM = 6*NS + 10
IF (ALLOCATED (INDX) ) DEALLOCATE (INDX)
ALLOCATE (INDX(3,NSM),STAT=IERR)
CALL AERR ('INDX(3,NSM)',IERR,INT(3*NSM))
if (jatri==3) then
call realloc(EDGEINDX, (/ 2,2*NSM /), keepExisting=.false., fill=0, stat=ierr)
call aerr('edgeindx(2,NSM)', ierr, int(2*NSM))
call realloc(TRIEDGE, (/ 3,NSM /), keepExisting=.false., fill=0, stat=ierr)
call aerr('triedge(3,NSM)', ierr, int(3*NSM))
else
call realloc(EDGEINDX, (/ 2,1 /), keepExisting=.false., fill=0, stat=ierr)
call realloc(TRIEDGE, (/ 3,1 /), keepExisting=.false., fill=0, stat=ierr)
end if
MAXTRI = NSM !?
numtri = NSM ! Input value should specify max nr of triangles in indx.
CALL TRICALL(jatri,XS,YS,NS,INDX,NUMTRI,EDGEINDX,NUMEDGE,TRIEDGE,XH,YH,NH,trisize)
if ( numtri.lt.0 ) nsm = -numtri
end do
call mess(LEVEL_INFO,'done')
END SUBROUTINE DLAUN
SUBROUTINE MAKEPLOTAREAS(NUMROW,NUMCOL,nsize)
implicit none
double precision :: dx
double precision :: dy
integer :: i, nsize
integer :: j
integer :: nsc
integer :: numcol
integer :: numrow
integer :: numsc
double precision :: x1sc
double precision :: x2sc
double precision :: xb
double precision :: xm
double precision :: xz
double precision :: y1sc
double precision :: y2sc
double precision :: yb
double precision :: ym
double precision :: yz
COMMON /GSCREENS/ X1SC(100),Y1SC(100),X2SC(100),Y2SC(100),NUMSC
NSC = 0
if (numrow == 1) then
yz = 0.4d0 ; yb = 0.8d0*(1d0-yz)
else
yz = 0.7d0 ; yb = 0.8d0*(1d0-yz)
endif
if (numcol < 3) then
xz = 0.7d0 ; xb = 0.5d0*(1d0-xz)
else
xz = 0.9d0 ; xb = 0.001d0 ! 05d0*(1d0-xz)
endif
if (nsize == 2) then
yz = 0.45d0 ; yb = 0.8d0*(1d0-yz)
endif
DY = yz / NUMROW
DX = xz / NUMCOL
XM = DX / 40
YM = DY / 40
DO 10 J = 1,NUMROW
DO 10 I = 1,NUMCOL
NSC = NSC + 1
X1SC(NSC) = xb + (I-1)*DX + XM
X2SC(NSC) = xb + (I )*DX - XM
Y1SC(NSC) = 1d0-yb - (J )*DY + YM
Y2SC(NSC) = 1d0-yb - (J-1)*DY - YM
10 CONTINUE
NUMSC = NSC
RETURN
END
!> bilineair interpolation between four nodes
subroutine bilin6(x, y, z, xp, yp, zp)
use m_missing
implicit none
double precision, dimension(2,2), intent(in) :: x, y !< node coordinates
double precision, dimension(2,2), intent(in) :: z !< node values
double precision, intent(in) :: xp, yp !< interpolant coordinates
double precision, intent(out) :: zp !< interpolant value
integer :: ierror
double precision :: A, B, C, D
double precision :: E, F, G, H
double precision :: P, Q, R, S
double precision :: xi1, eta1, xp1, yp1
double precision :: xi2, eta2, xp2, yp2
double precision :: aa, bb, cc, dis
double precision, parameter :: dtol = 1d-9
ierror = 1
zp = DMISS
! the interpolant has the form:
! xp(xi,eta) = A + B xi + C eta + D xi eta
! yp(xi,eta) = E + F xi + G eta + H xi eta
! zp(xi,eta) = P + Q xi + R eta + S xi eta
A = x(1,1)
B = x(2,1) - x(1,1) ! "xi-derivative"
C = x(1,2) - x(1,1) ! "eta-derivative"
D = (x(2,2) - x(2,1)) - (x(1,2)-x(1,1)) ! "xi-eta-derivative"
E = y(1,1)
F = y(2,1) - y(1,1) ! "xi-derivative"
G = y(1,2) - y(1,1) ! "eta-derivative"
H = (y(2,2) - y(2,1)) - (y(1,2)-y(1,1)) ! "xi-eta-derivative"
P = z(1,1)
Q = z(2,1) - z(1,1) ! "xi-derivative"
R = z(1,2) - z(1,1) ! "eta-derivative"
S = (z(2,2) - z(2,1)) - (z(1,2)-z(1,1)) ! "xi-eta-derivative"
! determine xi and eta from xp(xi,eta) = xp and yp(xi,eta)=yp
! aa xi^2 + bb xi + cc = 0
aa = D*F - B*H
bb = H*(xp-A) - B*G - D*(yp-E) + C*F
cc = G*(xp-A) - C*(yp-E)
dis = bb*bb - 4.0d0*aa*cc
if ( dis.lt.0d0 ) goto 1234
if ( abs(aa).gt.dtol ) then
dis = sqrt(dis)
aa = 0.5d0 / aa
xi1 = (-bb - dis) * aa
xi2 = (-bb + dis) * aa
else if ( abs(bb).gt.1d-9) then ! bb*xi + cc = 0
xi1 = -cc/bb
xi2 = xi1
else ! no or infinitely many solutions
goto 1234
end if
eta1 = C+D*xi1
if ( abs(eta1).gt.dtol ) then
eta1 = (xp-(A+B*xi1))/eta1
else
eta1 = G + H*xi2
if ( abs(eta1).gt.dtol ) then
eta1 = (yp-E-F*xi1)/ eta1
else
eta1 = DMISS
end if
end if
eta2 = C+D*xi2
if ( abs(eta2).gt.dtol ) then
eta2 = (xp-A-B*xi2)/eta2
else
eta2 = G + H*xi2
if ( abs(eta2).gt.dtol ) then
eta2 = (yp-E-F*xi2)/ eta2
else
eta2 = DMISS
end if
end if
! determine zp
do
! try first (xi,eta)
if ( xi1.ne.DMISS .and. eta1.ne.DMISS ) then
xp1 = A + B*xi1 + C*eta1 + D*xi1*eta1
yp1 = E + F*xi1 + G*eta1 + H*xi1*eta1
if ( abs(xp1-xp).lt.dtol .and. abs(yp1-yp).lt.dtol ) then
zp = P + Q*xi1 + R*eta1 + S*xi1*eta1
exit
end if
end if
! try second (xi,eta)
if ( xi2.ne.DMISS .and. eta2.ne.DMISS ) then
xp2 = A + B*xi2 + C*eta2 + D*xi2*eta2
yp2 = E + F*xi2 + G*eta2 + H*xi2*eta2
if ( abs(xp1-xp).lt.dtol .and. abs(yp1-yp).lt.dtol ) then
zp = P + Q*xi2 + R*eta2 + S*xi2*eta2
exit
end if
end if
! no solution found
goto 1234
end do
ierror = 0
! error handling
1234 continue
return
end
double precision function getdx(x1,y1,x2,y2)
use m_sferic
implicit none
double precision :: x1, y1, x2, y2
double precision :: xx1, yy1, xx2, yy2
double precision :: diff1, diff2
if (jsferic == 1) then
! fix for poles
diff1 = abs(abs(y1)-90d0)
diff2 = abs(abs(y2)-90d0)
if ( (diff1.le.dtol_pole .and. diff2.gt.dtol_pole) .or. &
(diff1.gt.dtol_pole .and. diff2.le.dtol_pole) ) then
getdx = 0d0
return
end if
xx1 = x1
xx2 = x2
if (xx1 - xx2 > 180d0) then
xx1 = xx1 - 360d0
else if (xx1 - xx2 < -180d0) then
xx1 = xx1 + 360d0
endif
xx1 = xx1*dg2rd
xx2 = xx2*dg2rd
yy1 = y1 *dg2rd
yy2 = y2 *dg2rd
csphi = dcos(0.5d0*(yy1+yy2))
getdx = ra*csphi*(xx2-xx1)
else
getdx = x2-x1
endif
end function getdx
double precision function getdy(x1,y1,x2,y2)
use m_sferic
implicit none
double precision :: x1, y1, x2, y2
double precision :: xx1, yy1,yy2
if (jsferic == 1) then
yy1 = y1*dg2rd
yy2 = y2*dg2rd
getdy = ra*(yy2-yy1)
else
getdy = y2-y1
endif
end function getdy
subroutine getdxdy(x1,y1,x2,y2,dx,dy)
use m_sferic
implicit none
double precision :: x1, y1, x2, y2, dx, dy, dx2, dy2, dum
double precision, external :: getdx, getdy
if (Jsferic == 1) then
if (jasferdistance == 1) then ! this is a fix
call sferdistance(x1,y1,x2,y1,dx)
if (x2 < x1) dx = -dx
call sferdistance(x1,y1,x1,y2,dy)
if (y2 < y1) dy = -dy
else
dx = getdx(x1,y1,x2,y2)
dy = getdy(x1,y1,x2,y2)
endif
else
dx = x2-x1
dy = y2-y1
endif
end subroutine getdxdy
subroutine sferdistance(x1,y1,x2,y2,d) ! == D3D's distance2
use m_sferic
implicit none
double precision :: x1,y1,x2,y2,d
double precision :: xx1, yy1, xx2, yy2, phi
xx1 = x1*dg2rd
xx2 = x2*dg2rd
yy1 = y1*dg2rd
yy2 = y2*dg2rd
phi = cos(yy1)*cos(yy2)*cos(xx1 - xx2) + sin(yy1)*sin(yy2)
d = ra*acos(phi)
end subroutine sferdistance
double precision function dprodin(x1,y1,x2,y2,x3,y3,x4,y4) ! inner product of two segments
use m_missing
implicit none
double precision :: x1,y1,x2,y2,x3,y3,x4,y4
double precision :: dx1,dy1,dx2,dy2
double precision :: getdx, getdy
dx1 = getdx(x1,y1,x2,y2)
dx2 = getdx(x3,y3,x4,y4)
dy1 = getdy(x1,y1,x2,y2)
dy2 = getdy(x3,y3,x4,y4)
dprodin = (dx1*dx2 + dy1*dy2)
return
end function dprodin
subroutine sincosdis(x1,y1,x2,y2,s,c,d) ! get sin, cos, length of a line segment
use m_missing
implicit none
double precision :: x1,y1,x2,y2,s,c,d
double precision :: dx1,dy1,dx2,dy2
double precision :: getdx, getdy
dx1 = getdx(x1,y1,x2,y2)
dy1 = getdy(x1,y1,x2,y2)
d = sqrt(dx1*dx1 + dy1*dy1)
if (d > 0d0) then
s = dy1/d
c = dx1/d
else
s = 0d0
c = 0d0
endif
end subroutine sincosdis
double precision function dprodout(x1,y1,x2,y2,x3,y3,x4,y4) ! out product of two segments
use m_missing
implicit none
double precision :: x1,y1,x2,y2,x3,y3,x4,y4
double precision :: dx1,dy1,dx2,dy2
double precision :: getdx, getdy
dx1 = getdx(x1,y1,x2,y2)
dx2 = getdx(x3,y3,x4,y4)
dy1 = getdy(x1,y1,x2,y2)
dy2 = getdy(x3,y3,x4,y4)
dprodout = (dx1*dy2 - dy1*dx2)
return
end function dprodout
!> Normalized inner product of two segments
!! NOTE that parallel lines may produce abs(dcosphi)=1+O(10^-16) > 1
!! in Debug builds, crashes subsequent acos calls! (not in Release)
double precision function dcosphi(x1,y1,x2,y2,x3,y3,x4,y4)
use m_missing
implicit none
double precision :: x1,y1,x2,y2,x3,y3,x4,y4
double precision :: dx1,dy1,dx2,dy2,r1,r2
double precision, external :: getdx, getdy
!call getdxdy(x1,y1,x2,y2,dx1,dy1)
!call getdxdy(x3,y3,x4,y4,dx2,dy2)
dx1 = getdx(x1,y1,x2,y2)
dx2 = getdx(x3,y3,x4,y4)
dy1 = getdy(x1,y1,x2,y2)
dy2 = getdy(x3,y3,x4,y4)
r1 = dx1*dx1+dy1*dy1
r2 = dx2*dx2+dy2*dy2
if (r1 == 0d0 .or. r2 == 0d0) then
dcosphi = dxymis
else
dcosphi = (dx1*dx2 + dy1*dy2)/sqrt(r1*r2)
endif
dcosphi = max( min(dcosphi,1d0) , -1d0)
return
end function dcosphi
double precision function dbdistance(x1,y1,x2,y2) ! distance point 1 -> 2
use m_missing
implicit none
double precision :: x1, y1, x2, y2
! locals
double precision :: ddx, ddy, rr
double precision :: getdx, getdy
if ( x1.eq.DMISS .or. x2.eq.DMISS .or. y1.eq.DMISS .or. y2.eq.DMISS ) then
dbdistance = 0d0
return
end if
ddx = getdx(x1,y1,x2,y2)
ddy = getdy(x1,y1,x2,y2)
rr = ddx*ddx + ddy*ddy
if (rr == 0d0) then
dbdistance = 0d0
else
dbdistance = sqrt(rr)
endif
end function dbdistance
SUBROUTINE dLINEDIS(X3,Y3,X1,Y1,X2,Y2,JA,DIS,XN,YN)
implicit none
integer :: ja
DOUBLE PRECISION :: X1,Y1,X2,Y2,X3,Y3,DIS,XN,YN
DOUBLE PRECISION :: R2,RL,X21,Y21,X31,Y31,getdx,getdy,dbdistance
! korste afstand tot lijnelement tussen eindpunten
JA = 0
X21 = getdx(x1,y1,x2,y2)
Y21 = getdy(x1,y1,x2,y2)
X31 = getdx(x1,y1,x3,y3)
Y31 = getdy(x1,y1,x3,y3)
R2 = dbdistance(x2,y2,x1,y1)
R2 = R2*R2
! IF (R2 .NE. 0) THEN
IF (R2 .GT. 1D-8) THEN
RL = (X31*X21 + Y31*Y21) / R2
RL = MAX( MIN(1d0,RL) , 0d0)
JA = 1
XN = X1 + RL*(x2-x1)
YN = Y1 + RL*(y2-y1)
DIS = dbdistance(x3,y3,xn,yn)
ELSE ! node 1 -> node 2
DIS = dbdistance(x3,y3,x1,y1)
ENDIF
RETURN
END subroutine DLINEDIS
SUBROUTINE dLINEDIS2(X3,Y3,X1,Y1,X2,Y2,JA,DIS,XN,YN,rl)
implicit none
integer :: ja
DOUBLE PRECISION :: X1,Y1,X2,Y2,X3,Y3,DIS,XN,YN
DOUBLE PRECISION :: R2,RL,X21,Y21,X31,Y31,getdx,getdy,dbdistance
! korste afstand tot lijnelement
JA = 0
X21 = getdx(x1,y1,x2,y2)
Y21 = getdy(x1,y1,x2,y2)
X31 = getdx(x1,y1,x3,y3)
Y31 = getdy(x1,y1,x3,y3)
R2 = dbdistance(x2,y2,x1,y1)
R2 = R2*R2
IF (R2 .NE. 0) THEN
RL = (X31*X21 + Y31*Y21) / R2
IF (0d0 .LE. RL .AND. RL .LE. 1d0) then
JA = 1
endif
XN = X1 + RL*(x2-x1)
YN = Y1 + RL*(y2-y1)
DIS = dbdistance(x3,y3,xn,yn)
ENDIF
RETURN
END subroutine DLINEDIS2
SUBROUTINE dLINEDIS3(X3,Y3,X1,Y1,X2,Y2,JA,DIS,XN,YN, RLOUT) ! 3: SORRY
implicit none
integer :: ja
DOUBLE PRECISION :: X1,Y1,X2,Y2,X3,Y3,DIS,XN,YN
DOUBLE PRECISION :: R2,RL,X21,Y21,X31,Y31,getdx,getdy,dbdistance
double precision :: RLout ! needed in orthogonalisenet/projection of boundary nodes
! korste afstand tot lijnelement tussen eindpunten
JA = 0
X21 = getdx(x1,y1,x2,y2)
Y21 = getdy(x1,y1,x2,y2)
X31 = getdx(x1,y1,x3,y3)
Y31 = getdy(x1,y1,x3,y3)
R2 = dbdistance(x2,y2,x1,y1)
R2 = R2*R2
RLout = 0d0
IF (R2 .NE. 0) THEN
RL = (X31*X21 + Y31*Y21) / R2
RLout = RL
RL = MAX( MIN(1d0,RL) , 0d0)
JA = 1
XN = X1 + RL*(x2-x1)
YN = Y1 + RL*(y2-y1)
DIS = dbdistance(x3,y3,xn,yn)
ENDIF
RETURN
END subroutine DLINEDIS3
!> Normalized vector in direction 1 -> 2
subroutine normalin(x1,y1,x2,y2,xn,yn)
use m_missing
implicit none
double precision :: x1, y1, x2, y2, xn, yn
! locals
double precision :: ddx, ddy, rr
double precision :: getdx, getdy
ddx = getdx(x1,y1,x2,y2)
ddy = getdy(x1,y1,x2,y2)
rr = ddx*ddx + ddy*ddy
if (rr == 0d0) then
xn = dxymis
yn = dxymis
else
rr = sqrt(rr)
xn = ddx / rr
yn = ddy / rr
endif
end subroutine normalin
!> Creates the relative unit normal vector to edge 1->2
!!
!! Vector is of unit length in Cartesian world.
!! Vector is almost unit length in spherical world, but its
!! x-component is scaled 1/cos(phi) such that in later uses:
!! (theta_B, theta_A) = (theta_A, phi_A) + alpha*(theta_n, phi_n)
!! the vectors 1->2 and A->B are perpendicular in Cartesian world,
!! not in spherical world. NOTE: the LENGTH of A->B in Cartesian
!! world implictly contains the earth radius and dg2rd already,
!! so make sure your alpha is in degrees.
subroutine normalout(x1,y1,x2,y2,xn,yn) ! normals out edge 1 2
use m_missing
use m_sferic
implicit none
double precision :: x1, y1, x2, y2, xn, yn
! locals
double precision :: ddx, ddy, rr
double precision :: getdx, getdy
ddx = getdx(x1,y1,x2,y2)
ddy = getdy(x1,y1,x2,y2)
rr = ddx*ddx + ddy*ddy
if (rr == 0d0) then
xn = dxymis
yn = dxymis
else
rr = sqrt(rr)
xn = ddy / rr
yn = -ddx / rr
endif
if (jsferic == 1) then
xn = xn / cos(dg2rd*0.5d0*(y1+y2) )
yn = yn
endif
end subroutine normalout
!> Computes the normal vector to a line 1-2, which is *outward* w.r.t.
!! an 'inside' point 3. Similar to normalout, except that the normal
!! vector may be flipped based on the 'inside' point.
subroutine normaloutchk(x1, y1, x2, y2, x3, y3, xn, yn, jaflip)
implicit none
double precision, intent(in) :: x1, y1 !< First point of line
double precision, intent(in) :: x2, y2 !< Second point of line
double precision, intent(in) :: x3, y3 !< Point that is considered 'inside'.
double precision, intent(out) :: xn, yn !< Output normal vector
integer, intent(out) :: jaflip !< Indicates whether normal was flipped (1) or not (0).
double precision :: din
double precision, external :: dprodin, dprodout
call normalout(x1,y1,x2,y2,xn,yn)
jaflip = 0
!din = dprodin(x1, y1, x1+xn, y1+yn, x1, y1, x3, y3)
!if (din > 0d0) then ! Check whether normal vector points really outward
! xn = -xn ! Using the previously stored internal point x4.
! yn = -yn
!end if
if ( dprodout(x1, y1, x1+xn, y1+yn, x1, y1, x2, y2)*dprodout(x1, y1, x3, y3, x1, y1, x2, y2) > 0d0 ) then
xn = -xn ! Using the previously stored internal point x4.
yn = -yn
jaflip = 1
endif
end subroutine normaloutchk
SUBROUTINE DUITPL(X1,Y1,X2,Y2,X3,Y3,X4,Y4,RU)
implicit none
double precision :: X1,Y1,X2,Y2,X3,Y3,X4,Y4,RU
double precision :: X12, y12, x34, y34
double precision :: GETDX, GETDY
X12 = GETDX(X1,Y1,X2,Y2)
Y12 = GETDY(X1,Y1,X2,Y2)
X34 = GETDX(X3,Y3,X4,Y4)
Y34 = GETDY(X3,Y3,X4,Y4)
RU = X12*Y34 - Y12*X34
RU = SIGN(1d0,RU)
RETURN
END SUBROUTINE DUITPL
!> Computes the enclosed area and length of a polygon.
!!
!! Only the first polygon is considered; whenever a missing value
!! is encountered, the polygon is 'closed'.
SUBROUTINE dAREAN( XX, YY, N, DAREA, DLENGTH, DLENMX )
USE M_MISSING
use m_sferic
implicit none
double precision, intent(in) :: XX(N), YY(N) !< Polygon points.
double precision, intent(out) :: DAREA !< Area enclosed within polygon.
double precision, intent(out) :: DLENGTH !< Length of polygon contour.
double precision, intent(out) :: DLENMX !< Length of longest segment in polygon contour.
integer, intent(in) :: n !< Nr. of polygon points.
integer :: i, iu, nend, jstart, jend
double precision :: DX, DY, Y0, DBDISTANCE, DLE, Y
DAREA = 0d0
DLENGTH = 0D0
Y0 = 1d30
NEND = 0
DLENMX = 0.D0
call get_startend(N,XX,YY,jstart,jend)
DO I = jstart,jend
IF (XX(I) .NE. dXYMIS) THEN
Y0 = MIN(Y0,YY(I))
NEND = I
ELSE
! dmiss encountered: end of first polygon.
! Only compute area for first polygon.
EXIT
ENDIF
ENDDO
DO I = jstart,jend
IU = I + 1
IF (IU .GT. NEND) IU = 1
IF (JSFERIC .EQ. 0) THEN
DX = ( XX(IU) - XX(I) )
Y = 0.5d0*(YY(IU)-Y0) + 0.5d0*(YY(I) -Y0)
ELSE
DX = ( XX(IU) - XX(I) )*DG2RD*RA*COS( DG2RD*( YY(IU)+YY(I) )/2 )
Y = RA*DG2RD*(0.5d0*( YY(IU) + YY(I) )-Y0)
ENDIF
DAREA = DAREA - DX*Y
DLE = DBDISTANCE( XX(I), YY(I), XX(IU), YY(IU) )
DLENGTH = DLENGTH + DLE
DLENMX = MAX (DLENMX, DLE)
ENDDO
DAREA = ABS(DAREA)
RETURN
END SUBROUTINE dAREAN
!> Write tecplot output in already opened file
!> note: file is closed when tim==tstop
subroutine tecplot_out(mtecfil, tim, Lwriheader)
use m_flowgeom
use m_flow
use m_flowtimes
use m_netw
implicit none
integer :: mtecfil
double precision, intent(in) :: tim
logical, intent(in) :: Lwriheader
integer :: i, j, k, icell
integer, dimension(4) :: inode
double precision, dimension(:), allocatable :: s1k, a1k, hk
character(len=256) :: zone
! Write header
if ( Lwriheader ) then
write(mtecfil, '(a)') 'Title = "output"'
write(mtecfil, '(a)') 'Variables = x y z e h u v'
end if
! Write cell centred data
! write (mtecfil, *) 'Zone f=point, t="', trim(zone),'"', &
! & ', i=', ndx
! do i = 1, ndx
! write (mtecfil, '(3f18.7)') xz(i), yz(i), s1(i) !, u1(inod)
! enddo
! Write net node data
! Interpolate s1 to net nodes
allocate(s1k(1:numk))
allocate(a1k(1:numk))
allocate(hk (1:numk))
s1k = 0.0d0
a1k = 0.0d0
hk = 0.0d0
! Loop over flow nodes (flow cells/elements) and over all its net nodes (corners) = expensive!
! Weighted interpolation based on cell area
do i = 1, ndxi
do j = 1, netcell(i)%N
k = netcell(i)%nod(j)
s1k(k) = s1k(k) + a1(i)*s1(i)
a1k(k) = a1k(k) + a1(i)
enddo
enddo
! Scale the interpolated water levels
s1k = s1k / a1k
! Compute total depth in net nodes
do i = 1, numk
hk(i) = max(0.0d0, s1k(i) - zk(i))
enddo
! Compute velocities to make sure the velocities in net nodes (corner points) have been computed
call setvelocityfield()
!
! Write zone information
!
write (zone , '(f16.4)' ) tim
write (mtecfil, '(3a)' ) 'Zone T ="', trim(zone),'"' !, numk
write (mtecfil, '(a,i0,a,i0,a)') 'N=', numk, ', E=', ndxi,', F=FEPOINT ET=QUADRILATERAL'
do i = 1, numk
write (mtecfil, '(7f18.7)') xk(i), yk(i), zk(i), s1k(i), hk(i), ucnx(i), ucny(i) !, u1k(i)
enddo
! Write connectivity table to file
do icell = 1, ndxi
do j=1,3
inode(j) = netcell(icell)%nod(j)
enddo
if (netcell(icell)%N == 3) then
inode(4) = netcell(icell)%nod(1)
else
inode(4) = netcell(icell)%nod(4)
endif
write (mtecfil, '(4(i0,3x))') inode(1), inode(2), inode(3), inode(4)
enddo
end subroutine tecplot_out
SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
implicit none
integer :: n,np,m,mp
double precision :: a,b
integer :: ipiv, indxr, indxc, i, j, k, L, LL, irow, icol
double precision :: big, dum, pivinv
! PARAMETER (NMAX=50)
! DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
DIMENSION A(NP,NP),B(NP,MP),IPIV(NP),INDXR(NP),INDXC(NP) ! SPvdP: set NMAX to N
DO 11 J=1,N
IPIV(J)=0
11 CONTINUE
DO 22 I=1,N
BIG=0.
DO 13 J=1,N
IF(IPIV(J).NE.1)THEN
DO 12 K=1,N
IF (IPIV(K).EQ.0) THEN
IF (ABS(A(J,K)).GE.BIG)THEN
BIG=ABS(A(J,K))
IROW=J
ICOL=K
ENDIF
ELSE IF (IPIV(K).GT.1) THEN
WRITE(*,*) 'Singular matrix'
ENDIF
12 CONTINUE
ENDIF
13 CONTINUE
IPIV(ICOL)=IPIV(ICOL)+1
IF (IROW.NE.ICOL) THEN
DO 14 L=1,N
DUM=A(IROW,L)
A(IROW,L)=A(ICOL,L)
A(ICOL,L)=DUM
14 CONTINUE
DO 15 L=1,M
DUM=B(IROW,L)
B(IROW,L)=B(ICOL,L)
B(ICOL,L)=DUM
15 CONTINUE
ENDIF
INDXR(I)=IROW
INDXC(I)=ICOL
IF (A(ICOL,ICOL).EQ.0.) WRITE(*,*) 'Singular matrix'
PIVINV=1./A(ICOL,ICOL)
A(ICOL,ICOL)=1.
DO 16 L=1,N
A(ICOL,L)=A(ICOL,L)*PIVINV
16 CONTINUE
DO 17 L=1,M
B(ICOL,L)=B(ICOL,L)*PIVINV
17 CONTINUE
DO 21 LL=1,N
IF(LL.NE.ICOL)THEN
DUM=A(LL,ICOL)
A(LL,ICOL)=0.
DO 18 L=1,N
A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18 CONTINUE
DO 19 L=1,M
B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19 CONTINUE
ENDIF
21 CONTINUE
22 CONTINUE
DO 24 L=N,1,-1
IF(INDXR(L).NE.INDXC(L))THEN
DO 23 K=1,N
DUM=A(K,INDXR(L))
A(K,INDXR(L))=A(K,INDXC(L))
A(K,INDXC(L))=DUM
23 CONTINUE
ENDIF
24 CONTINUE
RETURN
END
subroutine timdat(julday, timsec, idatum, itijd)
!!--copyright-------------------------------------------------------------------
! Copyright (c) 2006, WL | Delft Hydraulics. All rights reserved.
!!--disclaimer------------------------------------------------------------------
! This code is part of the Delft3D software system. WL|Delft Hydraulics has
! developed c.q. manufactured this code to its best ability and according to the
! state of the art. Nevertheless, there is no express or implied warranty as to
! this software whether tangible or intangible. In particular, there is no
! express or implied warranty as to the fitness for a particular purpose of this
! software, whether tangible or intangible. The intellectual property rights
! related to this software code remain with WL|Delft Hydraulics at all times.
! For details on the licensing agreement, we refer to the Delft3D software
! license and any modifications to this license, if applicable. These documents
! are available upon request.
!!--version information---------------------------------------------------------
! $Author$
! $Date$
! $Revision$
!!--description-----------------------------------------------------------------
!
! Function: returns date and time according actual time
! in minutes and itdate (julian day)
! in the form yyyymmdd hhmmss
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
!
implicit none
!
! Global variables
!
integer :: idatum ! Absolute date related to ITDATE and TIMSEC
integer , intent(out) :: itijd ! Absolute time related to ITDATE and TIMSEC
integer , intent(in) :: julday ! Description and declaration in inttim.igs
real(fp), intent(in) :: timsec ! Description and declaration in inttim.igs
!
! Local variables
!
integer(long) :: icurtm
integer :: iday
integer :: ihou
integer :: imin
integer :: imo
integer :: isec
integer :: iy
integer :: l
integer :: n
!
!! executable statements -------------------------------------------------------
!
! Define number of days, hours, minutes and seconds from TIMSEC
!
icurtm = nint(timsec,long)
iday = int(icurtm/86400_long)
icurtm = icurtm - int(iday,long)*86400_long
ihou = int(icurtm)/3600
icurtm = icurtm - int(ihou,long)*3600_long
imin = int(icurtm)/60
icurtm = icurtm - int(imin,long)*60_long
isec = int(icurtm)
!
! Convert julian day-number
!
l = julday + iday + 68569
n = 4*l/146097
l = l - (146097*n + 3)/4
iy = 4000*(l + 1)/1461001
l = l - 1461*iy/4 + 31
imo = 80*l/2447
iday = l - 2447*imo/80
l = imo/11
imo = imo + 2 - 12*l
iy = 100*(n - 49) + iy + l
!
! Define integer values for time and date
!
itijd = ihou*10000 + imin*100 + isec
idatum = abs(iy)*10000 + imo*100 + iday
idatum = sign(idatum, iy)
end subroutine timdat
!> Parses a manually preprocessed SVG file into 1D network 'drawing'.
!!
!! Format: each line should have one SVG command (m/M/c/l/z) with coordi nates.
subroutine parsekerst(filename)
use m_polygon
use m_missing
use unstruc_messages
use network_data
use unstruc_display, only:minmxns
use m_wearelt, only: rcir
implicit none
character(len=*), intent(in) :: filename
integer :: mfil, i, maxpts, it, maxt, KP, K1, LNU
character(len=2000) :: line
real, allocatable :: pts(:)
real :: startx, starty, curx, cury, x, y, t
double precision :: zp
maxpts = 200
allocate(pts(maxpts))
curx = 0.0
cury = 0.0
open(mfil, file=trim(filename))
do
read(mfil,'(a)',end=999) line
select case (line(1:1))
case ('c')
pts = dmiss
read(line(3:), *, end=97) pts
97 continue
do i = 1,maxpts-5,6
if (pts(i) == dmiss .or. pts(i+1) == dmiss) exit
!curx = curx + pts(i+4) ! Just take endpoint of bezier curve
!cury = cury + pts(i+5)
!npl = npl+1
!call increasepol(npl, 1)
!xpl(npl) = curx
!ypl(npl) = cury
maxt=5
do it=1,maxt
t = real(it)/real(maxt)
x = (1.0-t)**3 * (curx) &
+ 3.0*(1.0-t)**2 * t * (curx+pts(i)) &
+ 3.0*(1.0-t) * t**2 * (curx+pts(i+2)) &
+ t**3 * (curx+pts(i+4))
y = (1.0-t)**3 * (cury) &
+ 3.0*(1.0-t)**2 * t * (cury+pts(i+1)) &
+ 3.0*(1.0-t) * t**2 * (cury+pts(i+3)) &
+ t**3 * (cury+pts(i+5))
npl = npl+1
call increasepol(npl, 1)
xpl(npl) = x
ypl(npl) = y
end do
curx = x
cury = y
end do
case ('m','M') ! Move to new location
pts = dmiss
read(line(3:), *, end=98) pts
98 continue
if (line(1:1) == 'M') then
curx = pts(1)
cury = pts(2)
else
curx = curx + pts(1)
cury = cury + pts(2)
end if
startx = curx
starty = cury
npl = npl+1
call increasepol(npl+1, 1)
xpl(npl) = dmiss
ypl(npl) = dmiss
npl = npl+1
xpl(npl) = startx
ypl(npl) = starty
! If more than one point was given, treat as implicit subsequent lineto commands
do i = 3,maxpts-1,2
if (pts(i) == dmiss .or. pts(i+1) == dmiss) exit
npl = npl+1
call increasepol(npl, 1)
if (line(1:1) == 'M') then
curx = pts(i)
cury = pts(i+1)
else
curx = curx + pts(i)
cury = cury + pts(i+1)
end if
xpl(npl) = curx
ypl(npl) = cury
end do
case ('z') ! Close current subpath, return to subpath start
npl = npl+1
call increasepol(npl, 1)
xpl(npl) = startx
ypl(npl) = starty
curx = startx
cury = starty
case ('l')
pts = dmiss
read(line(3:), *, end=99) pts
99 continue
do i = 1,maxpts-1,2
if (pts(i) == dmiss .or. pts(i+1) == dmiss) exit
npl = npl+1
call increasepol(npl, 1)
curx = curx + pts(i)
cury = cury + pts(i+1)
xpl(npl) = curx
ypl(npl) = cury
end do
case default
! call mess('Unrecognized SVG command: '//trim(line))
end select
end do
999 continue
do i=1,npl
ypl(i) = -ypl(i)
end do
call zeronet()
KN3TYP=1
K1 = 0
rcir = 1d-10 ! isnode only for 'exact' hits
do i=1,npl
if (XPL(i) == dmiss) then
K1 = 0
cycle
end if
zp = 0d0
CALL ISNODE(KP, XPL(i), YPL(i), zp)
IF (KP .EQ. 0) THEN
CALL SETNEWPOINT(XPL(i), YPL(i), zp,KP)
ENDIF
IF (K1 .NE. 0) THEN
CALL CONNECTDBN(K1,KP,LNU)
END IF
K1=KP
end do
KN3TYP=2
call savepol() ! Store pol for optional restore
call delpol() ! Delete/hide the pol read from SVG file
CALL MINMXNS() ! New world scaling
!call wearel()
!CALL MERGENODESINPOLYGON() ! Merge 1D nodes that are too close.
deallocate(pts)
end subroutine parsekerst
subroutine getmdia(mdi) ! thanks herman
use unstruc_files
implicit none
integer :: mdi
mdi = mdia
end subroutine getmdia
subroutine setmdia(mdi) ! thanks herman, again
use unstruc_files
implicit none
integer :: mdi
mdia = mdi
end subroutine setmdia
!> make directory (also for linux)
subroutine makedir(dirname)
#ifdef __INTEL_COMPILER
use ifport
#endif
implicit none
character(len=*), intent(in) :: dirname
character(len=256) :: command
character(len=1), external :: get_dirsep
integer :: istat
logical :: l_exist
! write(6,"('Creating directory ', A128)") trim(dirname)
#ifdef __INTEL_COMPILER
inquire(directory = trim(dirname), exist = l_exist)
#else
! GNU
inquire(file = trim(dirname)//get_dirsep()//".", exist = l_exist)
#endif
if (l_exist) then
return
end if
if ( get_dirsep().eq.'/' ) then
! linux
command = "mkdir -p "//trim(dirname)
else
! windows
command = "mkdir "//trim(dirname)
! call iosDirMAKE(dirname)
end if
istat = system(command)
! Fortran2008, not available before Intel 15:
! call execute_command_line(command)
return
end subroutine
SUBROUTINE TIMLIN0()
implicit none
RETURN
END
SUBROUTINE FIRSTLIN(MRGF)
use unstruc_version_module, only: unstruc_version_full
implicit none
integer :: mrgf
CHARACTER TEX*80, RUNDAT*20
CALL DATUM(RUNDAT)
TEX = '* File creation date: ' //RUNDAT
WRITE(MRGF,'(A/A)') '* '//trim(unstruc_version_full), TEX
RETURN
END
!---------------------------------------------------------------
! the following subroutines use kdtree2
!---------------------------------------------------------------
!> find links crossed by polyline with kdtree2
subroutine find_crossed_links_kdtree2(treeinst,NPL,xpl,ypl,itype,nLinks,jaboundarylinks,numcrossedLinks, iLink, iPol, dSL, ierror)
!use m_polygon
use network_data, only: nump, numL, kn, xk, yk, lnn,lne
use m_flowgeom
use m_kdtree2
use m_sferic
use unstruc_messages
use m_missing
use m_alloc
implicit none
type(kdtree_instance), intent(inout) :: treeinst
integer, intent(in) :: NPL !< polyline length
double precision, dimension(NPL), intent(in) :: xpl, ypl !< polyline node coordinates
integer, intent(in) :: itype !< netlinks (1) or flowlinks(2)
integer, intent(in) :: nLinks !< number of links ( Lnx for flowlinks, numL for netlinks)
integer, intent(in) :: jaboundarylinks !< include boundary links (1) or not (0), flowlinks only
integer, intent(out) :: numcrossedLinks !< number of crossed flowlinks
integer, dimension(nLinks), intent(inout) :: iLink !< crossed flowlinks
integer, dimension(nLinks), intent(inout) :: iPol !< polygon section
double precision, dimension(nLinks), intent(inout) :: dSL !< polygon section cross location
integer, intent(out) :: ierror !< ierror (1) or not (0)
double precision, dimension(:), allocatable :: x, y
integer, dimension(:), allocatable :: ipolsection
double precision :: dmaxpollen, dmaxLinlen, dlinlen, R2search
integer :: num
integer, parameter :: jakdtree=1
integer, parameter :: MAXFIND=100
integer, parameter :: MINTREESIZE=0
double precision :: SL, SM, XCR, YCR, CRP
double precision :: xa, ya, xb, yb, af
integer :: i, k, L, N1, N2, NN, numnew
integer :: jacros, kint
integer :: LnxiORLnx
integer :: isactive
double precision, external :: dbdistance
ierror = 1
numcrossedLinks = 0
if ( NPL.lt.1 ) goto 1234 ! nothing to do
LnxiORLnx = 0
if ( itype.eq.1 ) then ! netlinks
LnxiORLnx = numL
else if ( itype.eq.2 ) then ! flowlinks
if ( jaboundarylinks.eq.1 ) then
LnxiORLnx = Lnx
else
LnxiORLnx = Lnxi
end if
end if
! allocate
allocate(ipolsection(NPL-1))
allocate(treeinst%qv(NTREEDIM))
! determine maximum polygon section length, and administer polygon sections
dmaxpollen = 0d0
num = 0
do i=1,NPL-1
if ( xpl(i).ne.DMISS .and. xpl(i+1).ne.DMISS ) then
num = num+1
ipolsection(num) = i
dmaxpollen = max(dmaxpollen, dbdistance(xpl(i),ypl(i),xpl(i+1),ypl(i+1)))
end if
end do
! check tree size and exit if the tree is too small
if ( num.lt.MINTREESIZE ) then
goto 1234
end if
! build kdtree
allocate(treeinst%sample_coords(NTREEDIM,NPL-1))
NN = min(MAXFIND,num)
! allocate
allocate(x(num), y(num))
! fill coordinates
do k=1,num
i=ipolsection(k)
x(k) = xpl(i)
y(k) = ypl(i)
end do
call build_kdtree(treeinst,num, x, y, ierror)
if ( ierror.ne.0 ) then
goto 1234
end if
! find crossed flowlinks
call mess(LEVEL_INFO, 'Finding crossed flowlinks...')
kint = max(LnxiORLnx/1000,1)
do L=1,LnxiORLnx
if (mod(L,kint) == 0) then
af = dble(L)/dble(LnxiORLnx)
call readyy('Finding crossed links', af)
endif
if ( itype.eq.1 ) then ! netlinks
call get_link_neighboringcellcoords(L,isactive,xa,ya,xb,yb)
if ( isactive.ne.1 ) then
cycle
end if
else if ( itype.eq.2 ) then ! flowlinks
n1 = ln(1,L) ; n2 = ln(2,L)
xa = xz(n1) ; ya = yz(n1)
xb = xz(n2) ; yb = yz(n2)
end if
! fill query vector
call make_queryvector_kdtree(treeinst,xa,ya)
! compute flowlink length
dlinlen = dbdistance(xa,ya,xb,yb)
! determine square search radius
R2search = 1.1d0*(dlinlen+dmaxpollen)**2 ! 1.1d0: safety
! count number of points in search area
NN = kdtree2_r_count(treeinst%tree,treeinst%qv,R2search)
if ( NN.eq.0 ) cycle ! no links found
! reallocate if necessary
call realloc_results_kdtree(treeinst,NN)
! find nearest NN points
call kdtree2_n_nearest(treeinst%tree,treeinst%qv,NN,treeinst%results)
jacros = 0
do i=1,NN
k = ipolsection(treeinst%results(i)%idx)
CALL CROSSinbox (XPL(k), YPL(k), XPL(k+1), YPL(k+1), Xa, Ya, Xb, Yb, jacros, SL, SM, XCR, YCR, CRP)
if ( jacros.eq.1 ) then
numcrossedLinks = numcrossedLinks + 1
if ( numcrossedLinks.gt.ubound(iLink,1) ) then
call mess(LEVEL_ERROR, 'find_crossed_links_kdtree2: array size too small')
end if
iLink(numcrossedLinks) = L
iPol(numcrossedLinks) = k
dSL(numcrossedLinks) = SL
end if
end do
end do
call readyy(' ', -1d0 )
call mess(LEVEL_INFO, 'done')
ierror = 0
1234 continue
! deallocate
if ( treeinst%itreestat.ne.ITREE_EMPTY ) call delete_kdtree2(treeinst)
if ( allocated(ipolsection) ) deallocate(ipolsection)
if ( allocated(x) ) deallocate(x)
if ( allocated(y) ) deallocate(y)
return
end subroutine find_crossed_links_kdtree2
!> build kdtree
subroutine build_kdtree(treeinst, N, x, y, ierror)
use m_kdtree2
use m_sferic
use unstruc_messages
use m_missing
use m_alloc
implicit none
type(kdtree_instance) :: treeinst
integer :: N !< number of entries
double precision, dimension(N) :: x, y !< coordinates
integer, intent(out) :: ierror !< error (1) or not (0)
integer :: k, num, ierr
ierror = 1
if ( N.eq.0 ) then
treeinst%itreestat = ITREE_EMPTY
goto 1234 ! nothing to do
end if
! build tree
call mess(LEVEL_INFO, 'Building kdtree...')
! check for spherical coordinates
if ( jsferic.eq.0 ) then
! Cartesian coordinates: 2D space
NTREEDIM = 2
else
! spherical coordinates: 3D space
NTREEDIM = 3
end if
if ( treeinst%itreestat.ne.ITREE_EMPTY ) call delete_kdtree2(treeinst)
allocate(treeinst%sample_coords(NTREEDIM,N),stat=ierr)
call aerr('sample_coords(NTREEDIM,N)', ierr, NTREEDIM*N)
treeinst%sample_coords = 0d0
! fill coordinates
num = 0
do k=1,N
if ( x(k).eq.DMISS .or. y(k).eq.DMISS ) then
cycle
end if
num = num+1
if ( jsferic.eq.0 ) then
! Cartesian coordinates: 2D space
treeinst%sample_coords(1,k) = x(k)
treeinst%sample_coords(2,k) = y(k)
else
! spherical coordinates: 3D space
treeinst%sample_coords(1,k) = Ra * cos(y(k)*dg2rd) * cos(x(k)*dg2rd)
treeinst%sample_coords(2,k) = Ra * cos(y(k)*dg2rd) * sin(x(k)*dg2rd)
treeinst%sample_coords(3,k) = Ra * sin(y(k)*dg2rd)
end if
end do
treeinst%tree => kdtree2_create(treeinst%sample_coords, rearrange=.true., sort=.true., dim=NTREEDIM)
! error handling
if ( kdtree2_ierror.ne.0 ) then
! call mess(LEVEL_DEBUG, 'kdtree error: kdtree2_ierror=', kdtree2_ierror)
call mess(LEVEL_info, 'kdtree error: kdtree2_ierror=', kdtree2_ierror)
goto 1234
end if
call mess(LEVEL_INFO, 'done')
! allocate query vector
allocate(treeinst%qv(NTREEDIM))
! allocate results array
IRESULTSIZE = INIRESULTSIZE
allocate(treeinst%results(IRESULTSIZE))
treeinst%itreestat = ITREE_OK
ierror = 0
1234 continue
return
end subroutine build_kdtree
!> make query vector for kdtree2
subroutine make_queryvector_kdtree(treeinst,x,y)
use m_sferic
use m_kdtree2
implicit none
double precision, intent(in) :: x, y !< coordinates
type(kdtree_instance) :: treeinst
! fill query vector
if ( jsferic.eq.0 ) then
treeinst%qv(1) = x
treeinst%qv(2) = y
else
treeinst%qv(1) = Ra * cos(y*dg2rd) * cos(x*dg2rd)
treeinst%qv(2) = Ra * cos(y*dg2rd) * sin(x*dg2rd)
treeinst%qv(3) = Ra * sin(y*dg2rd)
end if
return
end subroutine
!> resize results array
subroutine realloc_results_kdtree(treeinst,NN)
use m_kdtree2
implicit none
integer, intent(in) :: NN !< array size
type(kdtree_instance) :: treeinst
if ( .not.allocated(treeinst%results) ) then
IRESULTSIZE = INIRESULTSIZE
allocate(treeinst%results(IRESULTSIZE))
else
if ( NN.gt.IRESULTSIZE ) then
IRESULTSIZE = max(int(1.2d0*dble(NN))+1,INIRESULTSIZE)
deallocate(treeinst%results)
allocate(treeinst%results(IRESULTSIZE))
end if
end if
return
end subroutine
! delete kdtree
subroutine delete_kdtree2(treeinst)
use m_kdtree2
use unstruc_messages
implicit none
type(kdtree_instance) :: treeinst
if ( treeinst%itreestat.eq.ITREE_EMPTY ) return
call kdtree2_destroy(treeinst%tree)
if ( associated(treeinst%sample_coords) ) deallocate(treeinst%sample_coords)
if ( associated(treeinst%qv) ) deallocate(treeinst%qv)
if ( allocated(treeinst%results) ) deallocate(treeinst%results)
treeinst%itreestat = ITREE_EMPTY
IRESULTSIZE = 0
call mess(LEVEL_INFO, 'Delete kdtree...')
return
end subroutine delete_kdtree2
!> find nearest sample with kdtree2
subroutine find_nearest_sample_kdtree(treeinst,Ns,Ndim,xs,ys,zs,xk,yk,NN,isam,ierror)
use m_missing
use m_kdtree2
use m_alloc
use unstruc_messages
implicit none
integer, intent(in) :: Ns !< number of samples
integer, intent(in) :: Ndim !< dimension of sample vector
double precision, dimension(Ns), intent(in) :: xs, ys !< sample coordinates
double precision, dimension(Ndim,Ns), intent(in) :: zs !< sample values
double precision, intent(in) :: xk, yk !< query point coordinates
integer, intent(in) :: NN !< number of nearest samples
integer, dimension(NN), intent(out) :: isam !< nearest sample points
integer, intent(out) :: ierror !< error (>0), or not (0)
double precision, dimension(:), allocatable :: xx, yy
integer, parameter :: Nquerydim = 2 !< query vector length
integer :: i, num
type(kdtree_instance) :: treeinst
ierror = 0
! build tree if necessary
if ( treeinst%itreestat.ne.ITREE_OK ) then
if ( treeinst%itreestat.ne.ITREE_EMPTY ) then
call delete_kdtree2(treeinst)
end if
!! fill sample coordinates
! allocate(xx(NS), yy(NS))
! num = 0
! do i=1,NS
! if ( zs(1,i).ne.DMISS ) then ! check first dimension only
! num = num+1
! xx(num) = xs(i)
! yy(num) = ys(i)
! end if
! end do
!
! call build_kdtree(num,xx,yy,ierror)
! save memory by inserting samples directly (can give problems with DMISS-valued samples)
call build_kdtree(treeinst,NS,xs,ys,ierror)
!if ( allocated(xx) ) deallocate(xx)
!if ( allocated(yy) ) deallocate(yy)
end if
if ( ierror.eq.0 ) then
! fill query vector
call make_queryvector_kdtree(treeinst,xk,yk)
! reallocate if necessary
call realloc_results_kdtree(treeinst,NN)
! find nearest sample points
call kdtree2_n_nearest(treeinst%tree,treeinst%qv,NN,treeinst%results)
! copy to output
do i=1,NN
isam(i) = treeinst%results(i)%idx
end do
end if
return
end subroutine find_nearest_sample_kdtree
!> find flow cells with kdtree2
subroutine find_flowcells_kdtree(treeinst,Ns,xs,ys,inod,jaoutside,ierror)
use m_missing
use m_flowgeom
use m_kdtree2
use unstruc_messages
implicit none
type(kdtree_instance), intent(inout) :: treeinst
integer, intent(in) :: Ns !< number of samples
double precision, dimension(Ns), intent(in) :: xs, ys !< observation coordinates
double precision, dimension(:), allocatable :: xx, yy !< unique station coordinates
integer, dimension(:), allocatable :: iperm !< permutation array
integer, dimension(:), allocatable :: invperm !< inverse array
integer, dimension(Ns), intent(out) :: inod !< flow nodes
integer, intent(in) :: jaoutside !< allow outside cells (for 1D) (1) or not (0)
integer, intent(out) :: ierror !< error (>0), or not (0)
character(len=128) :: mesg
double precision :: dmaxsize, R2search, t0, t1
integer :: i, ip1, isam, in, k, N, NN
integer :: inum, num, jj
logical :: jadouble
double precision, external :: dbdistance
double precision :: dist_old, dist_new
ierror = 1
inod = 0
call klok(t0)
! reduce double stations ( see "fix for Wim" in subroutine rmdouble)
allocate(iperm(Ns), invperm(Ns))
iperm = 0
invperm = 0
allocate(xx(Ns),yy(Ns))
xx = 0d0
yy = 0d0
num = 0
do i=1,Ns
if ( xs(i).ne.DMISS .and. ys(i).ne.DMISS ) then
jadouble = .false.
do inum=1,num
if ( xs(i).eq.xx(inum) .and. ys(i).eq.yy(inum) ) then
jadouble = .true.
exit
end if
end do
if ( jadouble ) then
invperm(i) = -iperm(inum) ! store unique station index
cycle
end if
! new unique observation station
num = num+1
xx(num) = xs(i)
yy(num) = ys(i)
iperm(num) = i
invperm(i) = num ! really the inverse permutation
end if
end do
! build kdtree
call build_kdtree(treeinst, num, xx, yy, ierror)
if ( ierror.ne.0 ) then
goto 1234
end if
call mess(LEVEL_INFO, 'Finding flow nodes...')
! loop over flownodes
do k=1,Ndx
! fill query vector
call make_queryvector_kdtree(treeinst,xz(k),yz(k))
! compute maximum flowcell dimension
dmaxsize = 0d0
N = size(nd(k)%x)
do i=1,N
ip1=i+1; if ( ip1.gt.N ) ip1=ip1-N
dmaxsize = max(dmaxsize, dbdistance(nd(k)%x(i),nd(k)%y(i),nd(k)%x(ip1),nd(k)%y(ip1)))
end do
! determine square search radius
R2search = 1.1d0*dmaxsize**2 ! 1.1d0: safety
! count number of points in search area
NN = kdtree2_r_count(treeinst%tree,treeinst%qv,R2search)
if ( NN.eq.0 ) cycle ! no links found
! reallocate if necessary
call realloc_results_kdtree(treeinst,NN)
! find nearest NN samples
call kdtree2_n_nearest(treeinst%tree,treeinst%qv,NN,treeinst%results)
! check if samples are in cell
do i=1,NN
jj = treeinst%results(i)%idx
! find samples in original numbering (excluding the doubles)
isam = iperm(jj)
if ( k>ndx2D .and. kndx2D .and. k snap polygon to mesh
subroutine snappol(Nin, Xin, Yin, dsep, Nout, Xout, Yout, ipoLout, ierror)
use m_polygon
use m_missing
use m_alloc
use m_flowgeom
use network_data, only: xk, yk
implicit none
integer, intent(in) :: Nin !< thin-dyke polyline size
double precision, dimension(Nin), intent(in) :: Xin, Yin !< dsep-separated thin-dyke polyline coordinates
double precision, intent(in) :: dsep !< separator
integer, intent(out) :: Nout !< output polygon size
double precision, dimension(:), allocatable, intent(out) :: Xout, Yout !< output polygon coordinates, dim(Nout)
integer, dimension(:), allocatable, intent(out) :: ipoLout !< reference to input polyline (>0), seperator w.r.t. input polyline (0), dim(Nout)
integer, intent(out) :: ierror !< error (1) or not (0)
integer :: NumLinks, NDIM
double precision, dimension(:), allocatable :: dSL
integer, dimension(:), allocatable :: iLink, ipol
integer, dimension(:), allocatable :: ipolnr, indx
integer :: i, ii, iL, ipL, ipolsec, k1, k2, L, numpols
ierror = 1
Nout = 0
! save polygon
call savepol()
! allocate
allocate(iLink(Lnx))
iLink = 0
allocate(ipol(Lnx))
ipol = 0
allocate(dSL(Lnx))
dSL = 0d0
allocate(ipolnr(Nin))
ipolnr = 999
! number the input polyline segments
numpols = 0 ! polyline segment counter
i = 1 ! pointer in input array
do while ( i.lt.Nin )
! advance pointer to start of segment
do while ( (xin(i).eq.dsep .or. yin(i).eq.dsep) )
ipolnr(i) = 0
i=i+1
if ( i.gt.Nin ) exit
end do
! check for end of array
if ( i.gt.Nin ) exit
! mark this segment
numpols = numpols+1
do while ( xin(i).ne.dsep .and. yin(i).ne.dsep )
ipolnr(i) = numpols
i=i+1
if ( i.gt.Nin ) exit
end do
! check for end of array
if ( i.gt.Nin ) exit
end do
!! BEGIN DEBUG
! do i=1,NPL
! zpl(i) = dble(ipolnr(i))
! end do
! goto 1234
!! END DEBUG
! temporarily use output arrays to replace missing value
call realloc(xout, Nin, keepExisting=.false.)
call realloc(yout, Nin, keepExisting=.false.)
! replace missing values
do i=1,Nin
xout(i) = xin(i)
yout(i) = yin(i)
if ( xin(i).eq.dsep .or. yin(i).eq.dsep ) then
xout(i) = DMISS
yout(i) = DMISS
end if
end do
! snap polygon (note: xout and yout are temporary arrays)
call find_crossed_links_kdtree2(treeglob,Nin,xout,yout,2,Lnx,1,NumLinks,iLink, iPol, dSL, ierror)
if ( ierror.ne.0 .or. NumLinks.eq.0 ) goto 1234
! sort crossed flowlinks in increasing polyline order
allocate (indx(numLinks))
call indexxi(numLinks,iPol,indx)
! increase polygon array
call increasepol(3*NumLinks,0)
ii=1 ! pointer in indx array, sorted in increasing polygon number (iPol)
do ipL=1,numpols
! fill polygon with sections
i = 0
! advance pointer
do while ( ipolnr(iPol(indx(ii))).lt.ipL )
ii=ii+1
end do
if ( ii.gt.Numlinks ) then ! done
exit
end if
do while ( ipolnr(iPol(indx(ii))).eq.ipL )
!do iL=1,NumLinks
iL = indx(ii)
L = iLink(iL)
! check for matching polygon section
ipolsec = iPol(iL)
if ( ipolsec.lt.1 .or. ipolsec.ge.Nin ) then ! should not happen
continue
exit
end if
if ( ipolnr(ipolsec).eq.ipL .or. ipolnr(ipolsec+1).eq.ipL ) then
! find the netnodes
k1 = lncn(1,L)
k2 = lncn(2,L)
! fill the polyline array
i=i+1
xpl(i) = xk(k1)
ypl(i) = yk(k1)
i=i+1
xpl(i) = xk(k2)
ypl(i) = yk(k2)
i=i+1
xpl(i) = DMISS
ypl(i) = DMISS
else ! should not happen
continue
end if
ii = ii+1
if ( ii.gt.Numlinks ) exit ! done
end do
NPL = i
! merge polyline section parts
call merge_polylines()
if ( NPL.lt.2 ) cycle ! no polyline section found
! copy to output
NDIM = Nout + NPL + 1 ! add one for seperator
call realloc(Xout, NDIM, keepExisting=.true.)
call realloc(Yout, NDIM, keepExisting=.true.)
call realloc(ipoLout, NDIM, keepExisting=.true.)
do i=1,NPL
Xout(Nout+i) = xpl(i)
Yout(Nout+i) = ypl(i)
if ( xpl(i).eq.DMISS .or. xpl(i).eq.DMISS ) then
Xout(Nout+i) = dsep
Yout(Nout+i) = dsep
end if
ipoLout(Nout+i) = ipL
end do
Nout = Nout + NPL + 1 ! add one for seperator
Xout(Nout) = dsep
Yout(Nout) = dsep
ipoLout(Nout) = 0
if ( ii.gt.NumLinks ) exit ! done
end do
ierror = 0
1234 continue
if ( allocated(ipolnr) ) deallocate(ipolnr)
if ( allocated(iLink) ) deallocate(iLink)
if ( allocated(iPol) ) deallocate(iPol)
if ( allocated(dSL) ) deallocate(dSL)
call restorepol()
return
end subroutine snappol
!> snap point to flow node
subroutine snappnt(Nin, xin, yin, dsep, Nout, xout, yout, ipoLout, ierror)
use m_alloc
use m_flowgeom, only: xz, yz
implicit none
integer, intent(in) :: Nin !< thin-dyke polyline size
double precision, dimension(Nin) :: Xin, Yin !< dsep-separated thin-dyke polyline coordinates
double precision, intent(in) :: dsep !< missing value
integer, intent(out) :: Nout !< output polygon size
double precision, dimension(:), allocatable, intent(out) :: Xout, Yout !< output polygon coordinates, dim(Nout)
integer, dimension(:), allocatable, intent(out) :: ipoLout !< reference to input points (>0), no flownode found (0), dim(Nout)
integer, intent(out) :: ierror !< error (1) or not (0)
character(len=40), dimension(:), allocatable :: namobs
integer, dimension(:), allocatable :: kobs
integer :: i, k
integer :: jakdtree = 0
ierror = 1
Nout = 0
if ( Nin.lt.1 ) goto 1234
! allocate
allocate(namobs(Nin))
allocate(kobs(Nin))
do i=1,Nin
namobs(i) = ''
kobs(i) = 0
end do
call find_flownode(Nin, xin, yin, namobs, kobs, jakdtree, 1)
! copy to output
Nout = Nin
call realloc(xout, Nout, keepExisting=.false., fill=dsep)
call realloc(yout, Nout, keepExisting=.false., fill=dsep)
call realloc(ipoLout, Nout, keepExisting=.false., fill=0)
do i=1,Nout
k = kobs(i)
if ( k.gt.0 ) then
xout(i) = xz(k)
yout(i) = yz(k)
ipoLout(i) = i
else
xout(i) = dsep
yout(i) = dsep
ipoLout(i) = 0
end if
end do
ierror = 0
1234 continue
! deallocate
if ( allocated(namobs) ) deallocate(namobs)
if ( allocated(kobs) ) deallocate(kobs)
return
end subroutine snappnt
!> snap polyline to mesh boundary
!> 2D only
subroutine snapbnd(bndtype, Nin, Xin, Yin, dsep, Nout, Xout, Yout, ipoLout, ierror)
use timespace_triangle
use m_polygon
use m_missing
use network_data, only: kn, xk, yk, NumL, lne
use m_flowparameters, only: izbndpos
implicit none
character(len=*), intent(in) :: bndtype !< boundary condition type
integer, intent(in) :: Nin !< polyline size
double precision, dimension(Nin), intent(in) :: Xin, Yin !< dsep-separated polyline coordinates
double precision, intent(in) :: dsep !< separator
integer, intent(out) :: Nout !< output polygon size
double precision, dimension(:), allocatable, intent(out) :: Xout, Yout !< output polygon coordinates, dim(Nout)
integer, dimension(:), allocatable, intent(out) :: ipoLout !< reference to input polyline (>0), seperator w.r.t. input polyline (0), dim(Nout)
integer, intent(out) :: ierror !< error (1) or not (0)
double precision, dimension(:), allocatable :: xe, ye
double precision, dimension(:,:), allocatable :: xyen
double precision, dimension(:), allocatable :: xdum, ydum
integer, dimension(:), allocatable :: kce, ke, ki, kcs
integer, dimension(:), allocatable :: idx
double precision :: wL, wR
double precision :: xm, ym, crpm
double precision, dimension(4) :: xx, yy
double precision :: xzz, yzz, xci, yci, xce2, yce2
integer :: ierr, mx1Dend, Nx, numpols, jamiss
integer :: i, iend, j, k1, k2, k3, k4, kL, kR, L, m, num, NDIM
integer :: ja, isec, numsec, k, Lf
integer :: ioutput
integer, parameter :: INETLINKS = 0
integer, parameter :: IFLOWNODES = 1
integer, parameter :: IFLOWLINKS = 2
select case(trim(bndtype))
case( 'boundary' )
ioutput = INETLINKS
case( 'velocitybnd', 'dischargebnd' )
ioutput = IFLOWLINKS
case DEFAULT
ioutput = IFLOWNODES
end select
ierror = 1
Nout = 0
! save polygon
call savepol()
! count number of 2D links and 1D endpoints
call count_links(mx1Dend, Nx)
! allocate
allocate(xe(Nx), ye(Nx), xyen(2,Nx), kce(Nx), ke(Nin), ki(Nx))
allocate(kcs(Nin), xdum(Nin), ydum(Nin))
kce = 0
ke = 0
! replace missing values
do i=1,Nin
xdum(i) = xin(i)
ydum(i) = yin(i)
if ( xin(i).eq.dsep .or. yin(i).eq.dsep ) then
xdum(i) = DMISS
ydum(i) = DMISS
end if
end do
! make mirror cells (will set kce and ke)
call make_mirrorcells(Nx, xe, ye, xyen, kce, ke, ierror)
! set polyline mask
kcs = 1
! loop over the input polylines
numpols = 0 ! polyline counter
i = 1 ! pointer in input array
do while ( i.lt.Nin )
! advance pointer to start of segment
do while ( (xin(i).eq.dsep .or. yin(i).eq.dsep) )
i=i+1
if ( i.gt.Nin ) exit
end do
! check for end of array
if ( i.gt.Nin ) exit
! find end pointer of this polyline
numpols = numpols+1
iend = i
do while ( xin(iend).ne.dsep .and. yin(iend).ne.dsep )
iend=iend+1
if ( iend.gt.Nin ) exit
end do
iend = iend-1
!-------------------------------------------------------------------------------------------------------------------------
! the core of this subroutine (mostly copied from selectelset)
!-------------------------------------------------------------------------------------------------------------------------
! find boundary links
num = 0
do m = 1,Nx
if (iabs(kce(m)) == 1) then ! point is a possible candidate for a line boundary
call polyindexweight(xe(m), ye(m), xyen(1,m), xyen(2,m), Xdum(i:iend), Ydum(i:iend), kcs(i:iend), iend-i+1, kL, wL, kR, wR)
! if k1 > 0 this point can be dataprovided by this polyline
if (kL > 0 .or. kR > 0) then
if ( kce(m) .eq. -1 ) then
!errormessage = 'Boundary location already claimed; Overlap with other bnds?'
!return
continue
cycle
else
num = num + 1
ki(num) = m
kce(m) = -1 ! this tells you this point is already claimed by some bnd
endif
endif
endif
enddo
!-------------------------------------------------------------------------------------------------------------------------
! copy found net links to polygon segments
call increasepol(3*num, 0)
NPL = 0
do m=1,num
L = ki(m)
k1 = kn(1,L)
k2 = kn(2,L)
NPL = NPL+1
XPL(NPL) = xk(k1)
YPL(NPL) = yk(k1)
NPL = NPL+1
XPL(NPL) = xk(k2)
YPL(NPL) = yk(k2)
NPL = NPL+1
XPL(NPL) = DMISS
YPL(NPL) = DMISS
end do
! merge polyline section parts
call merge_polylines()
if ( ioutput.eq.IFLOWNODES .or. ioutput.eq.IFLOWLINKS ) then
! find flownodes or flowlinks for output
call realloc(idx, NPL, keepExisting=.false., fill=0)
numsec = 0
do i=1,num
m = ki(i)
! find polyline section (again)
call CROSSPOLY(xe(m),ye(m),xyen(1,m),xyen(2,m),XPL,YPL,NPL,xm,ym,crpm,ja,isec)
! remember which polyline segment points to this link
if ( isec.gt.0 ) then
if ( idx(isec).eq.0 ) then
numsec = numsec+1
idx(isec) = m
else ! should not happen
continue
end if
else ! should not happen
continue
end if
end do
end if
!-------------------------------------------------------------------------------------------------------------------------
! copy to output
NDIM = Nout + NPL + 1 ! add one for seperator
call realloc(Xout, NDIM, keepExisting=.true., fill = DMISS)
call realloc(Yout, NDIM, keepExisting=.true., fill = DMISS)
call realloc(ipoLout, NDIM, keepExisting=.true., fill = 0)
if ( ioutput.eq.INETLINKS ) then
do m=1,NPL
Xout(Nout+m) = xpl(m)
Yout(Nout+m) = ypl(m)
ipoLout(Nout+m) = numpols
end do
else if ( ioutput.eq.IFLOWLINKS ) then
num = 0 ! polyline size
jamiss = 0
do i=1,NPL
Xout(Nout+i) = DMISS
Yout(Nout+i) = DMISS
ipoLout(Nout+i) = numpols
L = idx(i)
if ( L.gt.0 .and. L.le.numL ) then
jamiss = 1
k1 = kn(1,L)
k2 = kn(2,L)
if ( k1.gt.0 .and. k2.gt.0 ) then
num = num+1
Xout(Nout+num) = 0.5d0*(xk(k1)+xk(k2))
Yout(Nout+num) = 0.5d0*(yk(k1)+yk(k2))
ipoLout(Nout+num) = numpols
end if
else
if ( jamiss.eq.1 ) then ! add seperator
num = num+1
Xout(Nout+num) = DMISS
Yout(Nout+num) = DMISS
ipolout(Nout+num) = numpols
jamiss = 0
end if
end if
end do
NPL = num
else if ( ioutput.eq.IFLOWNODES ) then
num = 0
jamiss = 0
do i=1,NPL
m = idx(i)
if ( m.gt.0 ) then
!-------------------------------------------------------------------------------------------------------------------------
! mostly copied from addexternalboundarypoints
!-------------------------------------------------------------------------------------------------------------------------
!if (kn(3,L) .ne. 1) then ! in 2D mirror cell
if ( m.le.numL ) then
jamiss = 1
L = m
k2 = iabs(lne(1,L))
k3 = kn(1,L); k4 = kn(2,L)
call mirrorcell(k2, xk(k3), yk(k3), xk(k4), yk(k4), xci, yci, xzz, yzz, xce2, yce2, xx, yy)
if (izbndpos == 0) then ! as in D3DFLOW
else if (izbndpos == 1) then ! on network boundary
xzz = 0.5d0*( xk(k3) + xk(k4 ) )
yzz = 0.5d0*( yk(k3) + yk(k4 ) )
else if (izbndpos == 2) then ! on specified boundary polyline
end if
num = num+1
Xout(Nout+num) = xzz
Yout(Nout+num) = yzz
ipoLout(Nout+num) = numpols
endif
!-------------------------------------------------------------------------------------------------------------------------
else
if ( jamiss.eq.1 ) then
num = num+1
Xout(Nout+num) = DMISS
Yout(Nout+num) = DMISS
ipoLout(Nout+num) = numpols
jamiss = 0
end if
end if
end do
NPL = num
end if
! remove trailing missing values
if ( NPL.gt.0 ) then
do while ( Xout(Nout+NPL).eq.DMISS .and. YOUT(Nout+NPL).eq.DMISS .and. NPL.gt.0)
NPL = NPL-1
if ( NPL.lt.1 ) exit
end do
end if
if ( NPL.gt.0 ) then
Nout = Nout + NPL + 1 ! add one for seperator
Xout(Nout) = DMISS
Yout(Nout) = DMISS
ipoLout(Nout) = 0
end if
! advance pointer
i = iend+1
! check for end of array
if ( i.gt.Nin ) exit
end do
! replace DMISS with seperator, if necessary
if ( dsep.ne.DMISS ) then
do i=1,Nout
if ( Xout(i).eq.DMISS .or. Yout(m).eq.DMISS ) then
Xout(i) = dsep
Yout(i) = dsep
end if
end do
end if
! trim output arrays to actual size
call realloc(Xout, Nout, keepExisting=.true.)
call realloc(Yout, Nout, keepExisting=.true.)
call realloc(ipolout, Nout, keepExisting=.true.)
ierror = 0
1234 continue
! deallocate
if ( allocated(xe) ) deallocate(xe)
if ( allocated(ye) ) deallocate(ye)
if ( allocated(xyen) ) deallocate(xyen)
if ( allocated(kce) ) deallocate(kce)
if ( allocated(ke) ) deallocate(ke)
if ( allocated(ki) ) deallocate(ki)
if ( allocated(kcs) ) deallocate(kcs)
if ( allocated(xdum) ) deallocate(xdum)
if ( allocated(ydum) ) deallocate(ydum)
if ( allocated(idx) ) deallocate(idx)
call restorepol()
return
end subroutine snapbnd
end module m_snappol
!< get jaopengl module variable
integer function iget_jaopengl()
use unstruc_opengl, only: jaopengl
implicit none
iget_jaopengl = jaopengl
return
end function iget_jaopengl
!< set jaopengl module variable
subroutine iset_jaopengl(jaopengl_loc)
use unstruc_opengl, only: jaopengl
use unstruc_model, only: md_jaopengl
implicit none
integer, intent(in) :: jaopengl_loc !< value to be set to jaopengl
if ( md_jaopenGL.eq.-1 ) then
jaopengl = jaopengl_loc
else
jaopengl = md_jaopengl
end if
return
end subroutine iset_jaopengl
!> prepares a matrix for solver test (as in "mpitest")
subroutine make_matrix(CFL, s1)
use m_reduce
use m_flowgeom
implicit none
double precision, intent(in) :: CFL !< CFL-number
double precision, dimension(Ndx), intent(in) :: s1 !< exact solution
double precision :: aufu
integer :: k1, k2, n, L
bbr = 1/CFL**2
ccr = 0d0
do L = 1,lnx
aufu = 1d0
k1 = ln(1,L)
k2 = ln(2,L)
bbr(k1) = bbr(k1) + aufu
bbr(k2) = bbr(k2) + aufu
ccr(Lv2(L)) = ccr(Lv2(L)) - aufu
enddo
! dd = 0d0
! do n=1,Ndx
! dd(n) = dd(n) + bb(n)*s1(n)
! end do
ddr = bbr*s1
do L=1,Lnx
k1 = ln(1,L)
k2 = ln(2,L)
ddr(k1) = ddr(k1) + ccr(Lv2(L))*s1(k2)
ddr(k2) = ddr(k2) + ccr(Lv2(L))*s1(k1)
enddo
return
end subroutine make_matrix
!> test iterative solver (as "mpitest")
subroutine soltest(iCFL,icgsolver_loc,maxsubmatvecs,iepsdiff,iepscg)
use m_partitioninfo
use m_timer
use unstruc_messages
use m_flowgeom
use network_data, only: xzw
use m_flowparameters
use m_reduce
use m_flow
use m_alloc
implicit none
integer, intent(in) :: iCFL !< wave-based Courant number
integer, intent(in) :: icgsolver_loc ! icgsolver (if > 0)
integer, intent(in) :: maxsubmatvecs ! maximum number of subiterations in Schwarz solver (if > 0)
integer, intent(in) :: iepsdiff ! -10log(tolerance in Schwarz iterations) (if > 0)
integer, intent(in) :: iepscg ! -10log(tolerance in inner iterations) (if > 0)
double precision, dimension(:), allocatable :: sex ! exact solution at cell centers
double precision :: CFL
double precision :: diffmax
integer :: NRUNS
integer :: i, ii, irun
integer :: ierror
CFL = 10d0
! maxdge = 0d0
! icgsolver = 4
! ipre = 0
Nruns = 1
if ( iCFL.gt.0d0 ) then
CFL = dble(iCFL)
end if
! settings from command line
if ( icgsolver_loc.gt.0 ) then
icgsolver = icgsolver_loc
end if
if ( iepsdiff.gt.0 ) then
epsdiff = 10d0**(-iepsdiff)
end if
if ( iepscg.gt.0 ) then
epscg = 10d0**(-iepscg)
end if
if ( maxmatvecs.gt.0 ) then
maxmatvecs = maxsubmatvecs
end if
call initimer()
call resetflow()
if ( jampi.eq.0 ) then
call flow_geominit(0)
if (Ndx == 0) then
call mess(LEVEL_INFO,'no network')
goto 1234
end if
else
call flow_geominit(1) ! first phase only
if ( Ndx.gt.0 ) then
if ( jatimer.eq.1 ) call starttimer(IPARTINIT)
call partition_init_2D('dummy', ierror) ! 2D only (hence the name, thanks to Herman for pointing this out)
if ( jatimer.eq.1 ) call stoptimer(IPARTINIT)
if ( ierror.ne.0 ) then
call mess(LEVEL_WARN,'Error in 2D partitioning initialization.')
goto 1234
end if
call update_geom(1) ! update geometry in ghost area
call flow_geominit(2) ! second phase
call update_geom(2) ! update geometry in ghost area
call disable_invalid_ghostcells_with_wu() ! disable ghost cells that are not being synchronised by setting wu's to zero
else
call mess(LEVEL_INFO,'no network')
goto 1234
end if
end if
call flow_allocflow()
! allocate solution
allocate(sex(Ndx))
! activate all cells
hu = epshu+1d0
! set exact solution
sex = xzw
if ( jatimer.eq.1 ) call starttimer(ITOTAL)
! prepare matrix
if ( jatimer.eq.1 ) call starttimer(IREDUCE)
call reducept(Ndx, Ndxi, Lnx)
if ( jatimer.eq.1 ) call stoptimer(IREDUCE)
! construct matrix and rhs
call make_matrix(CFL, sex)
! update overlapping ghost-parts of matrix
if ( jampi.eq.1 .and. jaoverlap.eq.1 ) then
call update_matrix(ierror)
end if
! pack matrix
call pack_matrix()
call realloc(ccrsav, ubound(ccr,1), lbound(ccr,1), keepExisting=.false., fill=0d0)
ccrsav = ccr
! solve system
do irun=1,Nruns
s1 = 0d0
ccr = ccrsav
! if (icgsolver.eq.6) call setPETSCmatrixEntries()
! call createPETSCPreconditioner(iprecond)
call solve_matrix(s1,Ndx,itsol)
end do
if ( jatimer.eq.1 ) call stoptimer(ITOTAL)
if ( jampi.eq.1 ) then
call update_ghosts(ITYPE_SALL,1,Ndx,s1,ierror)
end if
diffmax = 0d0
do i=1,Ndxi
if ( nd(i)%lnx.gt.0 ) then
if ( abs(s1(i)-sex(i)).gt.1d-10 ) then
continue
end if
diffmax = max(diffmax, abs(s1(i)-sex(i)))
end if
end do
do ii=1,nghostlist_sall(ndomains-1)
i = ighostlist_sall(ii)
if ( abs(s1(i)-sex(i)).gt.1d-10 ) then
continue
end if
end do
write(6,'("rank", I2, ", number of iterations: ", I4, ", max diff: ", E7.2)') my_rank, itsol, diffmax
if ( my_rank.eq.0 ) then
write(6,'(a,E8.2,a,E8.2)') ' WC-time solver [s]: ' , gettimer(1,ITOTALSOL), ' CPU-time solver [s]: ' , gettimer(0,ITOTALSOL)
write(6,'(a,E8.2,a,E8.2)') ' WC-time MPI comm [s]: ' , gettimer(1,IMPICOMM), ' CPU-time MPI comm [s]: ' , gettimer(0,IMPICOMM)
end if
! call mpi_barrier(DFM_COMM_DFMWORLD,ierr)
! call writemesg('Wallclock times')
! call printall(numt, t(3,:), tnams)
! call writemesg('CPU times')
! call printall(numt, tcpu(3,:), tnams)
1234 continue
if ( allocated(sex) ) deallocate(sex)
return
end subroutine soltest