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