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