subroutine SwanGradVel ( dep2, ux2, uy2, duxdx, duxdy, duydx, duydy, ivert ) ! ! --|-----------------------------------------------------------|-- ! | 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 ! ! 41.67: Marcel Zijlema ! ! Updates ! ! 41.67, August 2017: New subroutine ! ! Purpose ! ! Computes gradients of velocity components in vertex ! meant for computing transport velocities in spectral space ! ! Method ! ! Application of the Green-Gauss theorem with the assumption of ! a constant gradient over the controle volume (centroid dual) ! ! Modules used ! use ocpcomm4 use swcomm2 use swcomm3 use swcomm4 use SwanGriddata use SwanGridobjects ! implicit none ! ! Argument variables ! integer , intent(in) :: ivert ! counter corresponding to current vertex ! real, dimension(nverts), intent(in) :: dep2 ! water depth at current time level real, intent(out) :: duxdx ! derivative of ux2 to x real, intent(out) :: duxdy ! derivative of ux2 to y real, intent(out) :: duydx ! derivative of uy2 to x real, intent(out) :: duydy ! derivative of uy2 to y real, dimension(nverts), intent(in) :: ux2 ! ambient velocity in x-direction at current time level real, dimension(nverts), intent(in) :: uy2 ! ambient velocity in y-direction at current time level ! ! Local variables ! integer :: icell ! index of present cell integer, save :: ient = 0 ! number of entries in this subroutine integer :: jc ! loop counter integer :: jcell ! index of next cell ! integer, dimension(3) :: v ! vertices in present cell ! double precision :: area ! twices the area of centroid dual around present vertex real, dimension(MSC) :: arr ! auxiliary array real :: cslat ! cosine of latitude real :: dpmax ! maximum depth found in centroid dual real :: dpmin ! minimum depth found in centroid dual real, dimension(3) :: dloc ! local depth at vertices real :: dmaxc ! maximum depth found per cell of centroid dual real :: dminc ! minimum depth found per cell of centroid dual real, parameter :: drat= 5. ! ratio between maximum and minimum depths in centroid dual real :: ux0 ! u-velocity in centroid of present cell real :: ux1 ! u-velocity in centroid of next cell real :: uy0 ! v-velocity in centroid of present cell real :: uy1 ! v-velocity in centroid of next cell double precision :: x0 ! x-coordinate of the centroid of present cell double precision :: x1 ! x-coordinate of the centroid of next cell double precision :: y0 ! y-coordinate of the centroid of present cell double precision :: y1 ! y-coordinate of the centroid of next cell ! character(80) :: msgstr ! string to pass message ! 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,'SwanGradVel') ! duxdx = 0. duxdy = 0. duydx = 0. duydy = 0. ! ! if no frequency shift and no refraction, return ! if ( ITFRE == 0 .and. IREFR == 0 ) return ! ! point to vertex and cell objects ! vert => gridobject%vert_grid cell => gridobject%cell_grid ! if ( vert(ivert)%atti(VMARKER) == 1 ) return ! boundary vertex ! area = 0d0 dpmax = -99999. dpmin = 99999. ! ! loop over cells around considered vertex ! do jc = 1, vert(ivert)%noc ! ! get present 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) ! dloc(1) = dep2(v(1)) dloc(2) = dep2(v(2)) dloc(3) = dep2(v(3)) ! if ( dloc(1) <= DEPMIN .or. dloc(2) <= DEPMIN .or. dloc(3) <= DEPMIN ) goto 10 ! dminc = min ( min( dloc(1), dloc(2) ), dloc(3) ) if ( dminc < dpmin ) dpmin = dminc ! dmaxc = max ( max( dloc(1), dloc(2) ), dloc(3) ) if ( dmaxc > dpmax ) dpmax = dmaxc ! ! determine centroid of present cell ! x0 = cell(icell)%attr(CELLCX) y0 = cell(icell)%attr(CELLCY) ! ! determine u-velocity in centroid in present cell ! ux0 = ( ux2(v(1)) + ux2(v(2)) + ux2(v(3)) )/ 3. ! ! determine v-velocity in centroid in present cell ! uy0 = ( uy2(v(1)) + uy2(v(2)) + uy2(v(3)) )/ 3. ! ! get next cell in counterclockwise direction ! jcell = vert(ivert)%cell(jc)%atti(NEXTCELL) ! v(1) = cell(jcell)%atti(CELLV1) v(2) = cell(jcell)%atti(CELLV2) v(3) = cell(jcell)%atti(CELLV3) ! dloc(1) = dep2(v(1)) dloc(2) = dep2(v(2)) dloc(3) = dep2(v(3)) ! if ( dloc(1) <= DEPMIN .or. dloc(2) <= DEPMIN .or. dloc(3) <= DEPMIN ) goto 10 ! ! determine centroid of next cell ! x1 = cell(jcell)%attr(CELLCX) y1 = cell(jcell)%attr(CELLCY) ! ! determine u-velocity in centroid of next cell ! ux1 = ( ux2(v(1)) + ux2(v(2)) + ux2(v(3)) )/ 3. ! ! determine v-velocity in centroid in next cell ! uy1 = ( uy2(v(1)) + uy2(v(2)) + uy2(v(3)) )/ 3. ! ! compute contribution to area of centroid dual ! area = area + x0*y1 - x1*y0 ! ! compute contribution to x-gradient of velocity components ! duxdx = duxdx + ( ux0 + ux1 ) * real( y1 - y0 ) duydx = duydx + ( uy0 + uy1 ) * real( y1 - y0 ) ! ! compute contribution to y-gradient of velocity components ! duxdy = duxdy + ( ux0 + ux1 ) * real( x1 - x0 ) duydy = duydy + ( uy0 + uy1 ) * real( x1 - x0 ) ! enddo ! ! if ratio between max and min depth is too large, set gradients to zero and skip to next vertex ! !if ( dpmax > drat * dpmin ) goto 10 ! ! if area is non-positive, give error and go to next vertex ! if ( .not. area > 0d0 ) then write (msgstr, '(a,i5)') ' Area of centroid dual is negative or zero in vertex ', ivert call msgerr( 2, trim(msgstr) ) return endif ! duxdx = duxdx/real(area) duxdy = -duxdy/real(area) duydx = duydx/real(area) duydy = -duydy/real(area) ! ! in case of spherical coordinates, transform back to Cartesian coordinates ! if ( KSPHER > 0 ) then ! cslat = cos(DEGRAD*(vert(ivert)%attr(VERTY) + YOFFS)) ! duxdx = duxdx/(cslat * LENDEG) duxdy = duxdy/LENDEG duydx = duydx/(cslat * LENDEG) duydy = duydy/LENDEG ! endif ! return 10 duxdx = 0. duxdy = 0. duydx = 0. duydy = 0. ! end subroutine SwanGradVel