subroutine SwanInterpolateAc ( acintp, x, y, ac2, excpt )
!
!   --|-----------------------------------------------------------|--
!     | 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
!   40.90: Nico Booij
!
!   Updates
!
!   40.80, August 2007: New subroutine
!   40.90,   June 2008: improved interpolation near obstacles
!
!   Purpose
!
!   Interpolates action density to given point
!
!   Method
!
!   Look for closest vertex and determine triangle in which given point is located
!   Determine weighting coefficients for the corresponding vertices
!   Set weighting coeff to zero if there is an obstacle between given point and vertex
!   Interpolate action density using the resulting weighting coefficients
!
!   Modules used
!
    use ocpcomm4
    use swcomm2
    use swcomm3
    use m_obsta
    use SwanGriddata
    use SwanGridobjects
!
    implicit none
!
!   Argument variables
!
    real, dimension(MDC,MSC,nverts), intent(in) :: ac2    ! action density at current time
    real, dimension(MDC,MSC), intent(out)       :: acintp ! interpolated action density
    real, intent(in)                            :: x      ! x-coordinate of given point
    real, intent(in)                            :: y      ! y-coordinate of given point
    logical, intent(out)                        :: excpt  ! if true, value is undefined
!
!   Local variables
!
    integer                               :: icell     ! cell index
    integer, save                         :: ient = 0  ! number of entries in this subroutine
    integer, dimension(3)                 :: ivc       ! vertex indices in cyclic order
    integer                               :: ivert     ! vertex index
    integer                               :: jc        ! loop counter
    integer                               :: k         ! loop counter
    integer                               :: l         ! loop counter
    integer                               :: numcor    ! number of corner points in an obstacle
    integer, dimension(3)                 :: v         ! vertices in present cell
    !
    real                                  :: dxp       ! distance between given point and present vertex in x-direction
    real, dimension (3)                   :: dxv       ! difference of vertices of opposite side in x-coordinate
    real                                  :: dyp       ! distance between given point and present vertex in y-direction
    real, dimension (3)                   :: dyv       ! difference of vertices of opposite side in y-coordinate
    real                                  :: eps       ! a small number
    real                                  :: sumww     ! sum of the interpolation weights
    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                                  :: xb        ! user x-coordinate of begin of obstacle side
    real                                  :: xe        ! user x-coordinate of end of obstacle side
    real, dimension (3)                   :: xv        ! x-coordinate of the vertex
    real                                  :: yb        ! user y-coordinate of begin of obstacle side
    real                                  :: ye        ! user y-coordinate of end of obstacle side
    real, dimension (3)                   :: yv        ! y-coordinate of the vertex
    real, dimension (3)                   :: ww        ! weight of each vertex in the interpolation
    !
    character(80)                         :: msgstr    ! string to pass message
    !
    logical                               :: cellfound ! indicate whether cell containing given point is found or not
    logical, dimension (3)                :: cross     ! if true there is an obstacle between given point and vertex
    logical                               :: EQREAL    ! indicate whether two reals are equal or not
    logical                               :: obstcell  ! if true there is an obstacle in cell
    logical                               :: TCROSS    ! determines whether two line segments cross
    logical                               :: xonobst   ! not used
    !
    type(OBSTDAT), pointer                :: COBST     ! pointer to obstacle data
    !
    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,'SwanInterpolateAc')
    !
    ! point to vertex and cell objects
    !
    vert => gridobject%vert_grid
    cell => gridobject%cell_grid
    !
    ! initialize array for interpolated action density
    !
    acintp = 0.
    !
    excpt = .true.
    !
    ! find closest vertex for given point
    !
    call SwanFindPoint ( x, y, ivert )
    !
    ! if point not found, return
    !
    if ( ivert < 0 ) return
    !
    ! determine direction of given point to closest vertex
    !
    dxp = xcugrd(ivert) - x
    dyp = ycugrd(ivert) - y
    !
    ! if given point equals closest vertex, determine determine action density and return
    !
    if ( EQREAL(dxp,0.) .and. EQREAL(dyp,0.) ) then
       excpt  = .false.
       acintp = ac2(:,:,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
    !
    ! 2D linear interpolation on considered triangle is carried out now
    !
    excpt = .false.
    !
    !  get coordinates of the vertices
    !
    do k = 1, 3
       xv(k) = xcugrd(v(k))
       yv(k) = ycugrd(v(k))
       cross(k) = .false.
    enddo
    !
    ! determine difference in x and y of opposite side
    !
    do k = 1, 3
       ivc(2) = mod(k  ,3)+1
       ivc(3) = mod(k+1,3)+1
       dxv(k) = xv(ivc(3)) - xv(ivc(2))
       dyv(k) = yv(ivc(3)) - yv(ivc(2))
    enddo
    !
    ! determine whether there is an obstacle between given point and vertices
    !
    if ( NUMOBS > 0 ) then
       !
       COBST => FOBSTAC
       !
       do jc = 1, NUMOBS
          !
          numcor = COBST%NCRPTS
          if ( ITEST >= 120 ) write (PRINTF,10) jc, numcor
          !
          xb = COBST%XCRP(1)
          yb = COBST%YCRP(1)
          if ( ITEST >= 120 ) write (PRINTF,20) 1, xb+XOFFS, yb+YOFFS
          !
          do l = 2, numcor
             !
             xe = COBST%XCRP(l)
             ye = COBST%YCRP(l)
             if ( ITEST >= 120 ) write (PRINTF,20) l, xe+XOFFS, ye+YOFFS
             !
             ! loop over vertices
             !
             do k = 1, 3
                if ( TCROSS(x, xv(k), xb, xe, y, yv(k), yb, ye, xonobst) ) cross(k) = .true.
             enddo
             !
             xb = xe
             yb = ye
             !
          enddo
          !
          if (.not.associated(COBST%NEXTOBST)) exit
          COBST => COBST%NEXTOBST
          !
       enddo
       !
    endif
    !
    ! determine weighting coefficients
    !
    obstcell = .false.
    do k = 1, 3
       if (cross(k)) then
          ww(k) = 0.
          obstcell = .true.
       else
          ivc(1) = k
          ivc(2) = mod(k  ,3)+1
          ivc(3) = mod(k+1,3)+1
          ww(k) = ((x - xv(ivc(3))) * dyv(ivc(1)) - (y - yv(ivc(3))) * dxv(ivc(1))) / ( dxv(ivc(2)) * dyv(ivc(1)) - dyv(ivc(2)) * dxv(ivc(1)) )
       endif
    enddo
    if (obstcell) sumww = sum(ww)
    !
    ! use weighting coefficients to determine interpolated action density
    !
    do k = 1, 3
       if (obstcell) ww(k) = ww(k) / sumww
       acintp(:,:) = acintp(:,:) + ww(k) * ac2(:,:,v(k))
    enddo
    !
 10 format (' Obstacle number : ', i4,'  has ', i4, ' corners')
 20 format (' Corner number:', i4,'    Xp: ', e10.4, ' Yp: ', e11.4)
    !
end subroutine SwanInterpolateAc