subroutine SwanPropvelS ( cad , cas , ux2 , uy2 , & dep1 , dep2 , cax , cay , & kwave , cgo , spcsig, idcmin, & idcmax, ecos , esin , coscos, & sincos, sinsin, rdx , rdy , & mudl2 , jc ) ! ! --|-----------------------------------------------------------|-- ! | 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-2014 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.59: Erick Rogers ! 40.80: Marcel Zijlema ! 41.02: Marcel Zijlema ! 41.06: Gerbrant van Vledder ! 41.07: Marcel Zijlema ! 41.35: Casey Dietrich ! ! Updates ! ! 40.80, July 2007: New subroutine ! 40.59, August 2007: muddy bottom included ! 41.02, February 2009: adaption of velocities in case of diffraction ! 41.06, March 2009: add option of limitation of velocity in theta-direction ! 41.07, August 2009: add option of alternative formula for computation of ! wave transport velocity in theta-direction based on ! (x,y)-derivatives of the wave number ! (see Holthuijsen (2007), page 210, footnote 4) ! 41.35, March 2012: add option of limitation on csigma and ctheta ! ! Purpose ! ! computes wave transport velocities of energy in spectral space ! ! Modules used ! use ocpcomm4 use swcomm2 use swcomm3 use swcomm4 use m_diffr use SwanGriddata use SwanGridobjects use SwanCompdata ! implicit none ! ! Argument variables ! integer, intent(in) :: jc ! counter corresponding to current cell ! integer, dimension(MSC), intent(in) :: idcmax ! maximum frequency-dependent counter in directional space integer, dimension(MSC), intent(in) :: idcmin ! minimum frequency-dependent counter in directional space ! real, dimension(MDC,MSC), intent(out) :: cad ! wave transport velocity in theta-direction real, dimension(MDC,MSC), intent(out) :: cas ! wave transport velocity in sigma-direction real, dimension(MDC,MSC,ICMAX), intent(in) :: cax ! wave transport velocity in x-direction real, dimension(MDC,MSC,ICMAX), intent(in) :: cay ! wave transport velocity in y-direction real, dimension(MSC,ICMAX), intent(in) :: cgo ! group velocity real, dimension(MDC), intent(in) :: coscos ! help array containing cosine to power 2 of spectral directions real, dimension(nverts), intent(in) :: dep1 ! water depth at previous time level real, dimension(nverts), intent(in) :: dep2 ! water depth at current time level real, dimension(MDC), intent(in) :: ecos ! help array containing cosine of spectral directions real, dimension(MDC), intent(in) :: esin ! help array containing sine of spectral directions real, dimension(MSC,ICMAX), intent(in) :: kwave ! wave number real, dimension(nverts), intent(in) :: mudl2 ! mud layer at current time level real, dimension(MDC), intent(in) :: sincos ! help array containing sine * cosine of spectral directions real, dimension(MDC), intent(in) :: sinsin ! help array containing sine to power 2 of spectral directions 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 real, dimension(2), intent(in) :: rdx ! first component of contravariant base vector rdx(b) = a^(b)_1 real, dimension(2), intent(in) :: rdy ! second component of contravariant base vector rdy(b) = a^(b)_2 real, dimension(MSC), intent(in) :: spcsig ! relative frequency bins ! ! Local variables ! integer :: icell ! cell index integer :: id ! loop counter over direction bins integer :: iddum ! counter in directional space for considered sweep integer, save :: ient = 0 ! number of entries in this subroutine integer :: is ! loop counter over frequency bins integer :: iv1 ! first index in computational stencil integer :: iv2 ! second index in computational stencil integer :: iv3 ! third index in computational stencil integer :: jm ! counter corresponding to cell below considered sweep integer :: jp ! counter corresponding to cell above considered sweep integer :: k ! loop counter integer, dimension(3) :: v ! vertices in considered cell ! real :: alpha ! upper limit of CFL restriction real, dimension(3) :: cd ! coefficients for computing cad real, dimension(10) :: cs ! coefficients for computing cas real :: dhdx ! derivative of dep2 to x in considered sweep real :: dhdxm ! derivative of dep2 to x in sweep below considered sweep real :: dhdxp ! derivative of dep2 to x in sweep above considered sweep real :: dhdy ! derivative of dep2 to y in considered sweep real :: dhdym ! derivative of dep2 to y in sweep below considered sweep real :: dhdyp ! derivative of dep2 to y in sweep above considered sweep real, dimension(3) :: dloc ! local depth at vertices real, dimension(3) :: dm ! local mud layer at vertices real :: duxdx ! derivative of ux2 to x in considered sweep real :: duxdxm ! derivative of ux2 to x in sweep below considered sweep real :: duxdxp ! derivative of ux2 to x in sweep above considered sweep real :: duxdy ! derivative of ux2 to y in considered sweep real :: duxdym ! derivative of ux2 to y in sweep below considered sweep real :: duxdyp ! derivative of ux2 to y in sweep above considered sweep real :: duydx ! derivative of uy2 to x in considered sweep real :: duydxm ! derivative of uy2 to x in sweep below considered sweep real :: duydxp ! derivative of uy2 to x in sweep above considered sweep real :: duydy ! derivative of uy2 to y in considered sweep real :: duydym ! derivative of uy2 to y in sweep below considered sweep real :: duydyp ! derivative of uy2 to y in sweep above considered sweep real :: fac ! a factor real :: fac2 ! another factor real :: frlim ! frequency range in which limit on velocity in theta-direction is applied real :: kd ! help variable, wave number times water depth real :: pp ! power of the frequency dependent limiter on refraction real, dimension(2) :: rdxl ! first component of local contravariant base vector rdx(b) = a^(b)_1 real, dimension(2) :: rdyl ! second component of local contravariant base vector rdy(b) = a^(b)_2 ! real, dimension(MSC) :: arr ! auxiliary array real, dimension(MSC) :: dkdx ! derivative of wave number to x in considered sweep real, dimension(MSC) :: dkdxm ! derivative of wave number to x in sweep below considered sweep real, dimension(MSC) :: dkdxp ! derivative of wave number to x in sweep above considered sweep real, dimension(MSC) :: dkdy ! derivative of wave number to y in considered sweep real, dimension(MSC) :: dkdym ! derivative of wave number to y in sweep below considered sweep real, dimension(MSC) :: dkdyp ! derivative of wave number to y in sweep above considered sweep real, dimension(MSC,3) :: kloc ! local wave number ! 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,'SwanPropvelS') ! ! point to vertex and cell objects ! vert => gridobject%vert_grid cell => gridobject%cell_grid ! ! set velocities to zero ! cad = 0. cas = 0. ! ! if no frequency shift and no refraction, return ! if ( (ITFRE == 0 .or. (.not.DYNDEP .and. ICUR == 0)) .and. IREFR == 0 ) return ! ! initialize coefficients ! cd = 0. cs = 0. ! iv1 = vs(1) iv2 = vs(2) iv3 = vs(3) ! dloc(1) = dep2(iv1) dloc(2) = dep2(iv2) dloc(3) = dep2(iv3) ! ! if at least one vertex is dry, return ! if ( dloc(1) <= DEPMIN .or. dloc(2) <= DEPMIN .or. dloc(3) <= DEPMIN ) return ! ! compute the derivatives of the depth for the considered sweep ! if ( IREFR /= 0 .or. ICUR /= 0 ) then ! if ( IREFR == -1 ) then ! ! limitation of depths in upwave vertices ! dloc(2) = min( dep2(iv2), PNUMS(17)*dloc(1) ) dloc(3) = min( dep2(iv3), PNUMS(17)*dloc(1) ) ! endif ! dhdx = rdx(1) * (dloc(1)-dloc(2)) + rdx(2) * (dloc(1)-dloc(3)) dhdy = rdy(1) * (dloc(1)-dloc(2)) + rdy(2) * (dloc(1)-dloc(3)) ! endif ! ! compute the derivatives of the ambient current for the considered sweep ! if ( ICUR /= 0 ) then ! duxdx = rdx(1) * (ux2(iv1) - ux2(iv2)) + rdx(2) * (ux2(iv1) - ux2(iv3)) duxdy = rdy(1) * (ux2(iv1) - ux2(iv2)) + rdy(2) * (ux2(iv1) - ux2(iv3)) duydx = rdx(1) * (uy2(iv1) - uy2(iv2)) + rdx(2) * (uy2(iv1) - uy2(iv3)) duydy = rdy(1) * (uy2(iv1) - uy2(iv2)) + rdy(2) * (uy2(iv1) - uy2(iv3)) ! endif ! ! compute the derivatives of the wave number for the considered sweep ! if ( int(PNUMS(32)) == 1 ) then ! ! compute wave numbers for all frequencies ! call KSCIP1 (MSC, spcsig, dloc(1), kloc(1,1), arr, arr, arr) call KSCIP1 (MSC, spcsig, dloc(2), kloc(1,2), arr, arr, arr) call KSCIP1 (MSC, spcsig, dloc(3), kloc(1,3), arr, arr, arr) ! if ( IMUD == 1 ) then ! if (VARMUD) then dm(1) = mudl2(iv1) dm(2) = mudl2(iv2) dm(3) = mudl2(iv3) else dm(1) = PMUD(1) dm(2) = PMUD(1) dm(3) = PMUD(1) endif ! call KSCIP2 (MSC, spcsig, dloc(1), kloc(1,1), arr, arr, arr, arr, dm(1)) call KSCIP2 (MSC, spcsig, dloc(2), kloc(1,2), arr, arr, arr, arr, dm(2)) call KSCIP2 (MSC, spcsig, dloc(3), kloc(1,3), arr, arr, arr, arr, dm(3)) ! endif ! dkdx(:) = rdx(1) * (kloc(:,1)-kloc(:,2)) + rdx(2) * (kloc(:,1)-kloc(:,3)) dkdy(:) = rdy(1) * (kloc(:,1)-kloc(:,2)) + rdy(2) * (kloc(:,1)-kloc(:,3)) ! endif ! ! compute extra derivatives of depth, wave number and current belonging to neighbouring sweeps meant for refraction ! if ( IREFR /= 0 ) then ! if ( int(PNUMS(32)) == 0 ) then ! dhdxm = dhdx dhdym = dhdy dhdxp = dhdx dhdyp = dhdy ! else ! dkdxm = dkdx dkdym = dkdy dkdxp = dkdx dkdyp = dkdy ! endif ! if ( ICUR /= 0 ) then ! duxdxm = duxdx duxdym = duxdy duydxm = duydx duydym = duydy duxdxp = duxdx duxdyp = duxdy duydxp = duydx duydyp = duydy ! endif ! jm=vert(iv1)%noc-mod(vert(iv1)%noc+1-jc,vert(iv1)%noc) ! cell below considered sweep ! icell = vert(iv1)%cell(jm)%atti(CELLID) ! v(1) = cell(icell)%atti(CELLV1) v(2) = cell(icell)%atti(CELLV2) v(3) = cell(icell)%atti(CELLV3) ! ! compute local contravariant base vectors at present vertex in cell below considered sweep ! do k = 1, 3 if ( v(k) == iv1 ) then ! rdxl(1) = cell(icell)%geom(k)%rdx1 rdxl(2) = cell(icell)%geom(k)%rdx2 rdyl(1) = cell(icell)%geom(k)%rdy1 rdyl(2) = cell(icell)%geom(k)%rdy2 ! iv2 = v(mod(k ,3)+1) iv3 = v(mod(k+1,3)+1) ! exit endif enddo ! if ( KSPHER > 0 ) then do k = 1, 2 rdxl(k) = rdxl(k) / (COSLAT(1) * LENDEG) rdyl(k) = rdyl(k) / LENDEG enddo endif ! dloc(2) = dep2(iv2) dloc(3) = dep2(iv3) ! ! if at least one vertex is dry, skip and go to other cell ! if ( dloc(2) <= DEPMIN .or. dloc(3) <= DEPMIN ) goto 10 ! if ( IREFR == -1 ) then ! ! limitation of depths in upwave vertices ! dloc(2) = min( dep2(iv2), PNUMS(17)*dloc(1) ) dloc(3) = min( dep2(iv3), PNUMS(17)*dloc(1) ) ! endif ! if ( int(PNUMS(32)) == 0 ) then ! ! compute the derivatives of the depth in cell below considered sweep ! dhdxm = rdxl(1) * (dloc(1)-dloc(2)) + rdxl(2) * (dloc(1)-dloc(3)) dhdym = rdyl(1) * (dloc(1)-dloc(2)) + rdyl(2) * (dloc(1)-dloc(3)) ! else ! ! compute wave numbers for all frequencies ! call KSCIP1 (MSC, spcsig, dloc(2), kloc(1,2), arr, arr, arr) call KSCIP1 (MSC, spcsig, dloc(3), kloc(1,3), arr, arr, arr) ! if ( IMUD == 1 ) then ! if (VARMUD) then dm(2) = mudl2(iv2) dm(3) = mudl2(iv3) else dm(2) = PMUD(1) dm(3) = PMUD(1) endif ! call KSCIP2 (MSC, spcsig, dloc(2), kloc(1,2), arr, arr, arr, arr, dm(2)) call KSCIP2 (MSC, spcsig, dloc(3), kloc(1,3), arr, arr, arr, arr, dm(3)) ! endif ! ! compute the derivatives of the wave number in cell below considered sweep ! dkdxm(:) = rdxl(1) * (kloc(:,1)-kloc(:,2)) + rdxl(2) * (kloc(:,1)-kloc(:,3)) dkdym(:) = rdyl(1) * (kloc(:,1)-kloc(:,2)) + rdyl(2) * (kloc(:,1)-kloc(:,3)) ! endif ! ! compute the derivatives of the ambient current in cell below considered sweep ! if ( ICUR /= 0 ) then ! duxdxm = rdxl(1) * (ux2(iv1) - ux2(iv2)) + rdxl(2) * (ux2(iv1) - ux2(iv3)) duxdym = rdyl(1) * (ux2(iv1) - ux2(iv2)) + rdyl(2) * (ux2(iv1) - ux2(iv3)) duydxm = rdxl(1) * (uy2(iv1) - uy2(iv2)) + rdxl(2) * (uy2(iv1) - uy2(iv3)) duydym = rdyl(1) * (uy2(iv1) - uy2(iv2)) + rdyl(2) * (uy2(iv1) - uy2(iv3)) ! endif ! 10 continue ! jp = mod(jc,vert(iv1)%noc)+1 ! cell above considered sweep ! icell = vert(iv1)%cell(jp)%atti(CELLID) ! v(1) = cell(icell)%atti(CELLV1) v(2) = cell(icell)%atti(CELLV2) v(3) = cell(icell)%atti(CELLV3) ! ! compute local contravariant base vectors at present vertex in cell above considered sweep ! do k = 1, 3 if ( v(k) == iv1 ) then ! rdxl(1) = cell(icell)%geom(k)%rdx1 rdxl(2) = cell(icell)%geom(k)%rdx2 rdyl(1) = cell(icell)%geom(k)%rdy1 rdyl(2) = cell(icell)%geom(k)%rdy2 ! iv2 = v(mod(k ,3)+1) iv3 = v(mod(k+1,3)+1) ! exit endif enddo ! if ( KSPHER > 0 ) then do k = 1, 2 rdxl(k) = rdxl(k) / (COSLAT(1) * LENDEG) rdyl(k) = rdyl(k) / LENDEG enddo endif ! dloc(2) = dep2(iv2) dloc(3) = dep2(iv3) ! ! if at least one vertex is dry, skip ! if ( dloc(2) <= DEPMIN .or. dloc(3) <= DEPMIN ) goto 20 ! if ( IREFR == -1 ) then ! ! limitation of depths in upwave vertices ! dloc(2) = min( dep2(iv2), PNUMS(17)*dloc(1) ) dloc(3) = min( dep2(iv3), PNUMS(17)*dloc(1) ) ! endif ! if ( int(PNUMS(32)) == 0 ) then ! ! compute the derivatives of the depth in cell above considered sweep ! dhdxp = rdxl(1) * (dloc(1)-dloc(2)) + rdxl(2) * (dloc(1)-dloc(3)) dhdyp = rdyl(1) * (dloc(1)-dloc(2)) + rdyl(2) * (dloc(1)-dloc(3)) ! else ! ! compute wave numbers for all frequencies ! call KSCIP1 (MSC, spcsig, dloc(2), kloc(1,2), arr, arr, arr) call KSCIP1 (MSC, spcsig, dloc(3), kloc(1,3), arr, arr, arr) ! if ( IMUD == 1 ) then ! if (VARMUD) then dm(2) = mudl2(iv2) dm(3) = mudl2(iv3) else dm(2) = PMUD(1) dm(3) = PMUD(1) endif ! call KSCIP2 (MSC, spcsig, dloc(2), kloc(1,2), arr, arr, arr, arr, dm(2)) call KSCIP2 (MSC, spcsig, dloc(3), kloc(1,3), arr, arr, arr, arr, dm(3)) ! endif ! ! compute the derivatives of the wave number in cell above considered sweep ! dkdxp(:) = rdxl(1) * (kloc(:,1)-kloc(:,2)) + rdxl(2) * (kloc(:,1)-kloc(:,3)) dkdyp(:) = rdyl(1) * (kloc(:,1)-kloc(:,2)) + rdyl(2) * (kloc(:,1)-kloc(:,3)) ! endif ! ! compute the derivatives of the ambient current in cell above considered sweep ! if ( ICUR /= 0 ) then ! duxdxp = rdxl(1) * (ux2(iv1) - ux2(iv2)) + rdxl(2) * (ux2(iv1) - ux2(iv3)) duxdyp = rdyl(1) * (ux2(iv1) - ux2(iv2)) + rdyl(2) * (ux2(iv1) - ux2(iv3)) duydxp = rdxl(1) * (uy2(iv1) - uy2(iv2)) + rdxl(2) * (uy2(iv1) - uy2(iv3)) duydyp = rdyl(1) * (uy2(iv1) - uy2(iv2)) + rdyl(2) * (uy2(iv1) - uy2(iv3)) ! endif ! 20 continue ! endif ! ! compute wave transport velocity in sigma-direction ! if ( ITFRE /= 0 .and. (DYNDEP .or. ICUR /= 0) ) then ! ! compute time derivative of water depth ! if ( DYNDEP ) cs(2) = ( dep2(iv1) - dep1(iv1) ) * RDTIM ! ! compute coefficients depending on ambient currents ! if ( ICUR /= 0 ) then ! cs(3) = ux2(iv1) * dhdx cs(4) = uy2(iv1) * dhdy ! endif ! do is = 1, MSC ! ! compute frequency-dependent coefficients ! kd = min(30.,kwave(is,1) * dep2(iv1)) ! cs(1) = kwave(is,1) * spcsig(is) / sinh (2.* kd) cs(5) = -cgo(is,1) * kwave(is,1) if ( IDIFFR /= 0 ) cs(5) = cs(5)*DIFPARAM(iv1) ! cs( 6) = cs(1) * cs(2) cs( 7) = cs(1) * (cs(3)+cs(4)) cs( 8) = cs(5) * duxdx cs( 9) = cs(5) * (duxdy+duydx) cs(10) = cs(5) * duydy ! do iddum = idcmin(is), idcmax(is) id = mod ( iddum - 1 + MDC , MDC ) + 1 ! cas(id,is) = cs(6) if ( ICUR /= 0 ) cas(id,is) = cas(id,is) + cs(7) + coscos(id)*cs(8) + sincos(id)*cs(9) + sinsin(id)*cs(10) ! enddo enddo ! ! limit velocity using Courant number ! if ( int(PNUMS(33)) == 1 ) then ! alpha = PNUMS(34) ! do is = 1, MSC ! fac2 = alpha * FRINTF * spcsig(is) ! do iddum = idcmin(is), idcmax(is) id = mod ( iddum - 1 + MDC , MDC ) + 1 ! fac = fac2 * ( abs((rdx(1)+rdx(2))*cax(id,is,1)) + abs((rdy(1)+rdy(2))*cay(id,is,1)) ) ! if ( abs(cas(id,is)) > fac ) cas(id,is) = cas(id,is) * fac / abs(cas(id,is)) ! enddo ! enddo ! endif ! endif ! ! compute wave transport velocity in theta-direction ! if ( IREFR /= 0 ) then ! do is = 1, MSC ! ! compute frequency-dependent coefficients ! if ( int(PNUMS(32)) == 0 ) then ! kd = min(30.,kwave(is,1) * dep2(iv1)) ! cd(1) = spcsig(is) / sinh (2.* kd) ! cd(2) = cd(1) * dhdx cd(3) = cd(1) * dhdy ! else ! cd(1) = -cgo(is,1) / kwave(is,1) ! cd(2) = cd(1) * dkdx(is) cd(3) = cd(1) * dkdy(is) ! endif ! do iddum = idcmin(is), idcmax(is) id = mod ( iddum - 1 + MDC , MDC ) + 1 ! cad(id,is) = esin(id)*cd(2) - ecos(id)*cd(3) if ( IDIFFR /= 0 ) cad(id,is) = cad(id,is)*DIFPARAM(iv1) - DIFPARDX(iv1)*cgo(is,1)*esin(id) + DIFPARDY(iv1)*cgo(is,1)*ecos(id) if ( ICUR /= 0 ) cad(id,is) = cad(id,is) + sincos(id)*(duxdx-duydy) + sinsin(id)*duydx - coscos(id)*duxdy ! enddo ! if ( int(PNUMS(32)) == 0 ) then ! cd(2) = cd(1) * dhdxm cd(3) = cd(1) * dhdym ! else ! cd(2) = cd(1) * dkdxm(is) cd(3) = cd(1) * dkdym(is) ! endif ! id = mod ( idcmin(is) - 2 + MDC , MDC ) + 1 ! this direction belongs to sweep below considered sweep ! cad(id,is) = esin(id)*cd(2) - ecos(id)*cd(3) if ( IDIFFR /= 0 ) cad(id,is) = cad(id,is)*DIFPARAM(iv1) - DIFPARDX(iv1)*cgo(is,1)*esin(id) + DIFPARDY(iv1)*cgo(is,1)*ecos(id) if ( ICUR /= 0 ) cad(id,is) = cad(id,is) + sincos(id)*(duxdxm-duydym) + sinsin(id)*duydxm - coscos(id)*duxdym ! if ( int(PNUMS(32)) == 0 ) then ! cd(2) = cd(1) * dhdxp cd(3) = cd(1) * dhdyp ! else ! cd(2) = cd(1) * dkdxp(is) cd(3) = cd(1) * dkdyp(is) ! endif ! id = mod ( idcmax(is) + MDC , MDC ) + 1 ! this direction belongs to sweep above considered sweep ! cad(id,is) = esin(id)*cd(2) - ecos(id)*cd(3) if ( IDIFFR /= 0 ) cad(id,is) = cad(id,is)*DIFPARAM(iv1) - DIFPARDX(iv1)*cgo(is,1)*esin(id) + DIFPARDY(iv1)*cgo(is,1)*ecos(id) if ( ICUR /= 0 ) cad(id,is) = cad(id,is) + sincos(id)*(duxdxp-duydyp) + sinsin(id)*duydxp - coscos(id)*duxdyp ! enddo ! ! adapt velocity in case of spherical coordinates ! if ( KSPHER > 0 ) then ! cd(1) = tan(DEGRAD*(vert(iv1)%attr(VERTY) + YOFFS)) / REARTH ! do id = 1, MDC cd(2) = cd(1) * ecos(id) do is = 1, MSC cad(id,is) = cad(id,is) - cd(2)*(cax(id,is,1)*ecos(id) + cay(id,is,1)*esin(id)) enddo enddo ! endif ! ! limit velocity in some frequency range if requested ! if ( int(PNUMS(29)) == 1 ) then ! frlim = PI2*PNUMS(26) pp = PNUMS(27) ! do is = 1, MSC ! fac = min(1.,(spcsig(is)/frlim)**pp) ! do id = 1, MDC cad(id,is) = fac*cad(id,is) enddo enddo ! endif ! ! limit velocity using Courant number ! if ( int(PNUMS(35)) == 1 ) then ! alpha = PNUMS(36) ! fac2 = alpha * DDIR ! do is = 1, MSC ! do iddum = idcmin(is)-1, idcmax(is)+1 id = mod ( iddum - 1 + MDC , MDC ) + 1 ! fac = fac2 * ( abs((rdx(1)+rdx(2))*cax(id,is,1)) + abs((rdy(1)+rdy(2))*cay(id,is,1)) ) ! if ( abs(cad(id,is)) > fac ) cad(id,is) = cad(id,is) * fac / abs(cad(id,is)) ! enddo ! enddo ! endif ! endif ! end subroutine SwanPropvelS