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