module roelvink_module implicit none save !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! ! Jaap van Thiel de Vries, Robert McCall ! ! ! ! d.roelvink@unesco-ihe.org ! ! UNESCO-IHE Institute for Water Education ! ! P.O. Box 3015 ! ! 2601 DA Delft ! ! The Netherlands ! ! ! ! This library is free software; you can redistribute it and/or ! ! modify it under the terms of the GNU Lesser General Public ! ! License as published by the Free Software Foundation; either ! ! version 2.1 of the License, or (at your option) any later version. ! ! ! ! This library 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 ! ! Lesser General Public License for more details. ! ! ! ! You should have received a copy of the GNU Lesser General Public ! ! License along with this library; if not, write to the Free Software ! ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! ! USA ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! interface roelvink module procedure roelvink_1D module procedure roelvink_2D end interface roelvink interface baldock module procedure baldock_1D module procedure baldock_2D end interface baldock interface janssen_battjes module procedure janssen_battjes_1D module procedure janssen_battjes_2D end interface janssen_battjes contains subroutine roelvink_1D(par,s,km,i) use params use spaceparams use paramsconst implicit none type(spacepars) :: s type(parameters) :: par integer :: i,j real*8, dimension(s%ny+1) :: km,kmr,hr,H,arg ! Dissipation according to Roelvink (1993) H = sqrt(8.d0/par%rho/par%g*s%E(i,:)) hr = s%hh(i,:) + par%delta*s%H(i,:) kmr = min(max(km, 0.01d0), 100.d0) if (par%break /= BREAK_ROELVINK_DALY) then if (par%wci==1) then arg = -( H / (par%gamma*tanh(kmr*hr)/kmr))**par%n else arg = -( H / (par%gamma*hr ))**par%n endif s%Qb(i,:) = min(1.d0 - exp(max(arg,-100.d0)), 1.d0) else do j=1,s%ny+1 if (H(j) > par%gamma *hr(j)) s%Qb(i,j) = 1.d0 if (H(j) < par%gamma2*hr(j)) s%Qb(i,j) = 0.d0 enddo s%Qb(i,:) = max(s%Qb(i,:), 0.d0) endif s%D(i,:) = s%Qb(i,:) * 2.d0 * par%alpha * s%E(i,:) if (par%break == BREAK_ROELVINK1) then if (par%wci==1) then s%D(i,:) = s%D(i,:) * s%sigm(i,:)/2.d0/par%px; else s%D(i,:) = s%D(i,:) / par%Trep endif elseif (par%break == BREAK_ROELVINK2 .or. par%break == BREAK_ROELVINK_DALY) then ! Jaap: also wci switch for roelvink2 if (par%wci==1) then s%D(i,:) = s%D(i,:) * s%sigm(i,:)/2.d0/par%px * H/s%hh(i,:); else s%D(i,:) = s%D(i,:) / par%Trep * H/s%hh(i,:) endif end if end subroutine roelvink_1D subroutine roelvink_2D(par,s,km) use params use spaceparams implicit none type(spacepars) :: s type(parameters) :: par integer :: i real*8, dimension(s%nx+1,s%ny+1) :: km do i = 1,s%nx+1 call roelvink_1D(par,s,km(i,:),i) enddo end subroutine roelvink_2D subroutine baldock_1D(par,s,km,i) use params use spaceparams implicit none type(spacepars) :: s type(parameters) :: par integer :: i real*8, dimension(s%ny+1) :: km,kh,f,k,H,Hb,R,gamma ! Dissipation according to Baldock et al. (1998) if (par%wci==1) then f = s%sigm(i,:) / 2.d0 / par%px k = km else f = 1.d0 / par%Trep k = s%k(i,:) endif kh = s%k(i,:) * (s%hh(i,:) + par%delta*s%H(i,:)) if (par%wci == 1) then gamma = 0.76d0*kh + 0.29d0 !Jaap: spatial varying gamma according to Ruessink et al., 1998 else gamma = par%gamma endif H = sqrt(8.d0/par%rho/par%g*s%E(i,:)) Hb = tanh(gamma*kh/0.88d0)*(0.88d0/max(k,1d-10)) R = Hb/max(H,0.00001d0) s%Qb(i,:) = exp(-R**2) s%D (i,:) = 0.25d0 * par%alpha * f * par%rho * par%g * (Hb**2+H**2) * s%Qb(i,:) end subroutine baldock_1D subroutine baldock_2D(par,s,km) use params use spaceparams implicit none type(spacepars) :: s type(parameters) :: par integer :: i real*8, dimension(s%nx+1,s%ny+1) :: km do i = 1,s%nx+1 call baldock_1D(par,s,km(i,:),i) enddo end subroutine baldock_2D subroutine janssen_battjes_1D(par,s,km,i) use params use spaceparams use math_tools, only: xerf implicit none type(spacepars) :: s type(parameters) :: par integer :: i real*8 :: B real*8, dimension(s%ny+1) :: km,kh,f,k,H,Hb,R ! Dissipation according to Janssen and Battjes (2007) B = par%alpha if (par%wci==1) then f = s%sigm(i,:) / 2.d0 / par%px k = km else f = 1.d0 / par%Trep k = s%k(i,:) endif kh = s%k(i,:) * (s%hh(i,:) + par%delta*s%H(i,:)) H = sqrt(8.d0/par%rho/par%g*s%E(i,:)) Hb = tanh(par%gamma*kh/0.88d0)*(0.88d0/k) R = Hb/max(H,0.00001d0) s%Qb(i,:) = 1 + 4/(3*sqrt(par%px)) * (R**3 + 3/2*R) * exp(-R**2) - xerf(R) s%D (i,:) = 3*sqrt(par%px)/16 * B * f * par%rho * par%g * H**3/s%hh(i,:) * s%Qb(i,:) end subroutine janssen_battjes_1D subroutine janssen_battjes_2D(par,s,km) use params use spaceparams implicit none type(spacepars) :: s type(parameters) :: par integer :: i real*8, dimension(s%nx+1,s%ny+1) :: km do i = 1,s%nx+1 call janssen_battjes_1D(par,s,km(i,:),i) enddo end subroutine janssen_battjes_2D end module roelvink_module