real function SwanIntgratSpc ( p, fmin, fmax, spcsig, theta, wpar, ecs, uloc, vloc, acloc, itype )
!
! --|-----------------------------------------------------------|--
! | 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.87: Marcel Zijlema
!
! Updates
!
! 40.87, April 2008: New subroutine
!
! Purpose
!
! Determine p-th moment of energy density spectrum with respect to a part of frequency space
! (the integrand may contain power p of relative or absolute frequency, wave number or group velocity)
!
! Method
!
! Trapezoidal rule is applied
!
! Modules used
!
use ocpcomm4
use swcomm1
use swcomm3
!
implicit none
!
! Argument variables
!
integer, intent(in) :: itype ! indicate the type of integrand; represent power p of
! 1= relative frequency
! 2= absolute frequency
! 3= wave number
! 4= group velocity
!
real, intent(in) :: p ! power of the p-th moment
real, intent(in) :: fmax ! user-specified upper bound of frequency space for integration
real, intent(in) :: fmin ! user-specified lower bound of frequency space for integration
real :: uloc ! ambient u-velocity component at one location
real :: vloc ! ambient v-velocity component at one location
!
real, dimension(MDC,MSC), intent(in) :: acloc ! action density at one location
real, dimension(MDC), intent(in) :: ecs ! help array containing (co)sine of spectral directions
real, dimension(MSC), intent(in) :: spcsig ! relative frequency bins
real, dimension(MDC), intent(in) :: theta ! spectral directions
real, dimension(MSC), intent(in) :: wpar ! wave number or group velocity
!
! Local variables
!
integer :: id ! loop counter over direction bins
integer, save :: ient = 0 ! number of entries in this subroutine
integer :: is ! loop counter over frequency bins
!
real :: cia ! first coefficient in trapezoidal rule
real :: cib ! second coefficient in trapezoidal rule
real :: ctail ! coefficient for high frequency tail in case of (finite) upper bound
real :: ds ! width of frequency bin
real :: dsf ! remaining width of frequency bin near lower/upper bound
real :: ead ! auxiliary factor
real :: ftail ! factor for high frequency tail
real :: ehfr ! energy at high frequency
real :: omeg1 ! absolute frequency (upward bin)
real :: omeg2 ! another absolute frequency (current bin)
real :: pmom ! p-th moment of energy density spectrum
real :: q ! power of the p-th moment + 1
real :: us ! ambient velocity in wave direction
real :: wloc ! local wave parameter (wave number or group velocity)
!
! Structure
!
! Description of the pseudo code
!
! Source text
!
if (ltrace) call strace (ient,'SwanIntgratSpc')
!
q = p + 1.
!
pmom = 0.
!
if ( itype == 1 ) then
!
do is = 2, MSC
!
ds = spcsig(is) - spcsig(is-1)
!
if ( spcsig(is-1) >= fmin .and. spcsig(is) <= fmax ) then
!
cia = 0.5*spcsig(is-1)**q * ds
cib = 0.5*spcsig(is )**q * ds
!
elseif ( spcsig(is) > fmax ) then
!
dsf = fmax - spcsig(is-1)
cib = 0.5*fmax**q * dsf**2 / ds
cia = 0.5*(spcsig(is-1)**q + fmax**q)*dsf - cib
!
elseif ( spcsig(is) > fmin ) then
!
dsf = spcsig(is) - fmin
cia = 0.5*fmin**q * dsf**2 / ds
cib = 0.5*(spcsig(is)**q + fmin**q)*dsf - cia
!
else
!
cia = 0.
cib = 0.
!
endif
!
do id = 1, MDC
!
ead = cia * acloc(id,is-1) + cib * acloc(id,is)
pmom = pmom + ead * ecs(id) * DDIR
!
enddo
!
if ( spcsig(is) > fmax ) exit
!
enddo
!
! add tail contribution, if appropriate
!
if ( fmax > spcsig(MSC) ) then
!
if ( MSC > 3 ) then
!
if ( PWTAIL(1) <= q ) then
call msgerr (2, ' power of moment is too large compared to tail power')
SwanIntgratSpc = -999.
return
endif
!
ftail = 1. / (PWTAIL(1) - q)
!
if ( fmax > 100. ) then
!
ctail = 0.
!
else
!
ctail = fmax**q * (spcsig(MSC)/fmax)**PWTAIL(1)
!
endif
!
do id = 1, MDC
!
ehfr = acloc(id,MSC) * spcsig(MSC) * ecs(id)
pmom = pmom + ehfr * (spcsig(MSC)**q - ctail) * DDIR * ftail
!
enddo
!
endif
!
endif
!
elseif ( itype == 2 ) then
!
do id = 1, MDC
!
us = uloc*cos(theta(id)+ALCQ) + vloc*sin(theta(id)+ALCQ)
!
do is = 2, MSC
!
ds = spcsig(is) - spcsig(is-1)
!
omeg1 = spcsig(is-1) + wpar(is-1) * us
omeg2 = spcsig(is ) + wpar(is ) * us
!
if ( spcsig(is-1) >= fmin .and. spcsig(is) <= fmax ) then
!
cia = 0.5*omeg1**p * spcsig(is-1) * ds
cib = 0.5*omeg2**p * spcsig(is ) * ds
!
elseif ( spcsig(is) > fmax ) then
!
dsf = fmax - spcsig(is-1)
omeg2 = fmax + (wpar(is)*dsf+wpar(is-1)*(ds-dsf)) * us / ds
cib = 0.5*omeg2**p * fmax * dsf**2 / ds
cia = 0.5*(omeg1**p * spcsig(is-1) + omeg2**p * fmax)*dsf - cib
!
elseif ( spcsig(is) > fmin ) then
!
dsf = spcsig(is) - fmin
omeg1 = fmin + (wpar(is-1)*dsf+wpar(is)*(ds-dsf)) * us / ds
cia = 0.5*omeg1**p * fmin * dsf**2 / ds
cib = 0.5*(omeg2**p * spcsig(is) + omeg1**p * fmin)*dsf - cia
!
else
!
cia = 0.
cib = 0.
!
endif
!
ead = cia * acloc(id,is-1) + cib * acloc(id,is)
pmom = pmom + ead * ecs(id) * DDIR
!
if ( spcsig(is) > fmax ) exit
!
enddo
!
enddo
!
! add tail contribution, if appropriate
!
if ( fmax > spcsig(MSC) ) then
!
if ( MSC > 3 ) then
!
ftail = 1. / (PWTAIL(1) - 1.)
!
if ( fmax > 100. ) then
!
ctail = 0.
!
else
!
ctail = fmax**q * (spcsig(MSC)/fmax)**PWTAIL(1)
!
endif
!
do id = 1, MDC
!
us = uloc*cos(theta(id)+ALCQ) + vloc*sin(theta(id)+ALCQ)
omeg2 = spcsig(MSC) + wpar(MSC) * us
!
ehfr = acloc(id,MSC) * spcsig(MSC) * ecs(id)
pmom = pmom + ehfr * (omeg2**p * spcsig(MSC) - ctail) * DDIR * ftail
!
enddo
!
endif
!
endif
!
elseif ( itype == 3 .or. itype == 4 ) then
!
do is = 2, MSC
!
ds = spcsig(is) - spcsig(is-1)
!
if ( spcsig(is-1) >= fmin .and. spcsig(is) <= fmax ) then
!
cia = 0.5*wpar(is-1)**p * spcsig(is-1) * ds
cib = 0.5*wpar(is )**p * spcsig(is ) * ds
!
elseif ( spcsig(is) > fmax ) then
!
dsf = fmax - spcsig(is-1)
wloc = (wpar(is)*dsf+wpar(is-1)*(ds-dsf)) / ds
cib = 0.5*wloc**p * fmax * dsf**2 / ds
cia = 0.5*(wpar(is-1)**p * spcsig(is-1) + wloc**p * fmax)*dsf - cib
!
elseif ( spcsig(is) > fmin ) then
!
dsf = spcsig(is) - fmin
wloc = (wpar(is-1)*dsf+wpar(is)*(ds-dsf)) / ds
cia = 0.5*wloc**p * fmin * dsf**2 / ds
cib = 0.5*(wpar(is)**p * spcsig(is) + wloc**p * fmin)*dsf - cia
!
else
!
cia = 0.
cib = 0.
!
endif
!
do id = 1, MDC
!
ead = cia * acloc(id,is-1) + cib * acloc(id,is)
pmom = pmom + ead * ecs(id) * DDIR
!
enddo
!
if ( spcsig(is) > fmax ) exit
!
enddo
!
! add tail contribution, if appropriate
!
if ( fmax > spcsig(MSC) ) then
!
if ( MSC > 3 ) then
!
ftail = 1. / (PWTAIL(1) - 1.)
!
if ( fmax > 100. ) then
!
ctail = 0.
!
else
!
if ( itype == 3 ) then
!
ctail = (1./GRAV)**p * fmax**(2*p+1) * (spcsig(MSC)/fmax)**PWTAIL(1)
!
elseif ( itype == 4 ) then
!
ctail = (0.5*GRAV)**p * fmax**(1-p) * (spcsig(MSC)/fmax)**PWTAIL(1)
!
endif
!
endif
!
do id = 1, MDC
!
ehfr = acloc(id,MSC) * spcsig(MSC) * ecs(id)
pmom = pmom + ehfr * (wpar(MSC)**p * spcsig(MSC) - ctail) * DDIR * ftail
!
enddo
!
endif
!
endif
!
endif
!
SwanIntgratSpc = pmom
!
end function SwanIntgratSpc