subroutine SwanInterpolatePoint ( foutp, x, y, finp, excval )
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! Authors
!
! 40.80: Marcel Zijlema
!
! Updates
!
! 40.80, August 2007: New subroutine
!
! Purpose
!
! Interpolates given scalar to given point
!
! Method
!
! First, look for closest vertex and next, interpolate given scalar inside triangle where given point is resided
!
! Modules used
!
use ocpcomm4
use swcomm2
use swcomm3
use SwanGriddata
use SwanGridobjects
!
implicit none
!
! Argument variables
!
real, intent(in) :: excval ! exception value for given scalar
real, dimension(nverts), intent(in) :: finp ! given scalars defined on the computational grid
real, intent(out) :: foutp ! output scalar at given point
real, intent(in) :: x ! x-coordinate of given point
real, intent(in) :: y ! y-coordinate of given point
!
! Local variables
!
integer :: icell ! cell index
integer, save :: ient = 0 ! number of entries in this subroutine
integer :: ivert ! vertex index
integer :: jc ! loop counter
integer :: k ! loop counter
integer, dimension(3) :: v ! vertices in present cell
!
real :: carea ! area of the present cell
real :: dxp ! distance between given point and present vertex in x-direction
real :: dyp ! distance between given point and present vertex in y-direction
real :: eps ! a small number
real :: phi1 ! value of given scalar in first vertex of considered cell
real :: phi2 ! value of given scalar in second vertex of considered cell
real :: phi3 ! value of given scalar in third vertex of considered cell
real :: phic ! value of given scalar in centroid of considered cell
real :: th ! direction of given point to present vertex
real :: th1 ! direction of one face pointing to present vertex
real :: th2 ! direction of another face pointing to present vertex
real :: thdiff ! difference between th and th2
real, dimension(2) :: vec12 ! translation vector of coordinates: vertex2 - vertex1
real, dimension(2) :: vec23 ! translation vector of coordinates: vertex3 - vertex2
real, dimension(2) :: vec31 ! translation vector of coordinates: vertex1 - vertex3
real :: xc ! x-coordinate of the cell-centroid
real :: yc ! y-coordinate of the cell-centroid
real :: xgrs ! x-component of gradient scalar vector
real :: ygrs ! y-component of gradient scalar vector
!
character(80) :: msgstr ! string to pass message
!
logical :: cellfound ! indicate whether cell containing given point is found or not
logical :: EQREAL ! indicate whether two reals are equal or not
!
type(celltype), dimension(:), pointer :: cell ! datastructure for cells with their attributes
type(verttype), dimension(:), pointer :: vert ! datastructure for vertices with their attributes
!
! Structure
!
! Description of the pseudo code
!
! Source text
!
if (ltrace) call strace (ient,'SwanInterpolatePoint')
!
! point to vertex and cell objects
!
vert => gridobject%vert_grid
cell => gridobject%cell_grid
!
! assign exception value to output scalar (possibly overwritten by interpolated value)
!
foutp = excval
!
! find closest vertex for given point
!
call SwanFindPoint ( x, y, ivert )
!
! if point not found, give warning and return
!
if ( ivert < 0 ) then
write (msgstr, '(a,f12.4,a,f12.4,a)') ' Point (',x+XOFFS,',',y+YOFFS,') not given in computational grid'
call msgerr( 1, trim(msgstr) )
return
endif
!
! if exception value found in closest vertex, return
!
if ( EQREAL(finp(ivert),excval) ) return
!
! determine direction of given point to closest vertex
!
dxp = xcugrd(ivert) - x
dyp = ycugrd(ivert) - y
!
! if given point equals closest vertex, determine output quantity and return
!
if ( EQREAL(dxp,0.) .and. EQREAL(dyp,0.) ) then
foutp = finp(ivert)
return
endif
!
th = atan2(dyp,dxp)
!
cellfound = .false.
!
! loop over cells around closest vertex
!
celloop: do jc = 1, vert(ivert)%noc
!
! get cell and its vertices
!
icell = vert(ivert)%cell(jc)%atti(CELLID)
!
v(1) = cell(icell)%atti(CELLV1)
v(2) = cell(icell)%atti(CELLV2)
v(3) = cell(icell)%atti(CELLV3)
!
! get directions of faces to closest vertex
!
do k = 1, 3
if ( v(k) == ivert ) then
th1 = cell(icell)%geom(k)%th1
th2 = cell(icell)%geom(k)%th2
exit
endif
enddo
!
thdiff = th - th2
do
if ( abs(thdiff) <= PI ) exit
th = th - sign (2., thdiff) * PI
thdiff = th - th2
enddo
!
! is given point inside considered cell?
!
if ( vert(ivert)%atti(VMARKER) == 1 ) then ! boundary vertex
eps = PI/360.
else
eps = 0.
endif
!
if ( th > th1-eps .and. th <= th2+eps ) then
cellfound = .true.
exit celloop
endif
!
enddo celloop
!
! if cell containing given point not found, give warning and return
!
if ( .not.cellfound ) then
write (msgstr, '(a,f12.4,a,f12.4,a)') ' No triangle containing point (',x+XOFFS,',',y+YOFFS,') is found'
call msgerr( 1, trim(msgstr) )
return
endif
!
! determine output scalar in vertices
!
phi1 = finp(v(1))
phi2 = finp(v(2))
phi3 = finp(v(3))
!
! 2D linear interpolation on considered triangle is carried out only if there are no exception values
!
if ( .not.EQREAL(phi1,excval) .and. .not.EQREAL(phi2,excval) .and. .not.EQREAL(phi3,excval) ) then
!
! determine centroid and area of found cell
!
xc = cell(icell)%attr(CELLCX )
yc = cell(icell)%attr(CELLCY )
carea = cell(icell)%attr(CELLAREA)
!
! determine output scalar in centroid
!
phic = ( phi1 + phi2 + phi3 ) / 3.
!
! determine translation vectors of found cell
!
vec12(1) = xcugrd(v(2)) - xcugrd(v(1))
vec12(2) = ycugrd(v(2)) - ycugrd(v(1))
vec23(1) = xcugrd(v(3)) - xcugrd(v(2))
vec23(2) = ycugrd(v(3)) - ycugrd(v(2))
vec31(1) = xcugrd(v(1)) - xcugrd(v(3))
vec31(2) = ycugrd(v(1)) - ycugrd(v(3))
!
! determine gradient scalar vector inside found cell based on outward normals
! Note: the outward normal is obtained by rotating the translation vector
! over 90 degrees in clockwise direction
!
xgrs = vec23(2)*phi1 + vec31(2)*phi2 + vec12(2)*phi3
ygrs = -vec23(1)*phi1 - vec31(1)*phi2 - vec12(1)*phi3
!
xgrs = -0.5*xgrs/carea
ygrs = -0.5*ygrs/carea
!
! determine output scalar inside considered triangle by means of 2D interpolation
! using constant gradient scalar vector
!
foutp = phic + xgrs*(x - xc) + ygrs*(y - yc)
!
endif
!
end subroutine SwanInterpolatePoint