module sandmud_module
use params
! public routines
!
!public readsandmud
implicit none
public sandmud_erodep
public erosed
public eromud
public erosand
public sand_mud
public vanRijn84
public shld
public initsed
public sandmud_sb_vr
public tau_bed
public settling
contains
!===================================================================================!
subroutine sandmud_erodep(s,par)
!
use params
use spaceparams
use bedcomposition_module
use readbedcomposition_module, only: bc
use message_module
use handlearray_module
!
implicit none
!
include 'sedparams.inc'
!
type(parameters) :: par
type(spacepars) :: s
type(message_stack), pointer :: messages ! message stack
!
real*8, dimension((par%nx+1)*(par%ny+1),par%ngd) :: sm_ero ! [kg/m2/s] erosion rate
real*8, dimension((par%nx+1)*(par%ny+1),par%ngd) :: sm_depo_ex ! [m/s] volumetric deposition flux per unit surface
real*8, dimension((par%nx+1)*(par%ny+1),par%ngd) :: sm_Csedeq ! [kg/m3] equilibrium concentration
real*8, dimension(par%ngd,par%nx+1,par%ny+1) :: ws ! [m/s] fall velocity
real*8, dimension(par%nx+1,par%ny+1) :: tau_tot ! [N/m2] bottom shear stress
real*8, dimension(par%nx+1,par%ny+1) :: U ! [m/s] velocity magnitude at cell center
real*8, dimension(par%nx+1,par%ny+1) :: Chezy ! [m/1/2)/s] Chezy coefficient
real*8, dimension(par%ngd,(par%nx+1)*(par%ny+1)) :: sink ! [m/s] deposition rate
real*8, dimension(par%ngd,(par%nx+1)*(par%ny+1)) :: sinkf ! [m/s] deposition rate to fluff layer
real*8, dimension(par%ngd,(par%nx+1)*(par%ny+1)) :: sour ! [kg/m2/s] erosion flux => [m3/m2/s] before exit the subroutine
real*8, dimension(par%ngd,(par%nx+1)*(par%ny+1)) :: sourf ! [kg/m2/s] erosion flux from fluff layer
real*8, dimension(par%ngd,(par%nx+1)*(par%ny+1)) :: rsedeqb ! [m3/m3] equilibrium bed concentration
real*8 :: Te,kvis,Sster,cc1,cc2,wster,w
integer :: i,j,jg,pos
!
!
bc%settings%nmlb = 1
bc%settings%nmub = (par%nx+1)*(par%ny+1)
!
! Initialize variables
s%ero = 0.d0
s%depo_ex = 0.d0
s%ceqbg = 0.d0
!
! Compute fall velocity with expression from van Rijn (date?)
! TODO: find a version for cohesive sediment
call settling(s,par,ws,bc%settings%rhofrac)
!
! Compute velocity magnitude at each cell
U = s%ue**2+s%ve**2
!
! Compute bottom shear stress at each cell
! TODO: add effect of turbulence due to wave breaking ?
call tau_bed(s,par,tau_tot,s%hh,U,s%urms)
!
! Provide Chezy coefficient to all cells
Chezy = par%C
!
! Determine erosion and deposition fluxes
call erosed(bc , bc%settings%nmlb , bc%settings%nmub, par%dt , par%morfac , &
& ws , U , s%hh , Chezy , tau_tot , &
& bc%settings%nfrac, bc%settings%rhofrac, par%D50 , par%D90 , bc%settings%sedtyp, &
& sink , sinkf , sour , sourf , messages , &
& rsedeqb , s%R , s%cg , s%kb , s%kturb , &
& s%tbore , s%cf , s%wetz , par%sedcal , s%urms , &
& par )
!
!
sm_ero = transpose(sour)
sm_depo_ex = transpose(sink)
sm_Csedeq = transpose(rsedeqb) ! Equilibrium bed concentration (from bedload)
do jg=1,bc%settings%nfrac
sm_Csedeq(:,jg) = sm_Csedeq(:,jg)/bc%settings%rhofrac(jg)
sm_ero(:,jg) = sm_ero(:,jg)/bc%settings%rhofrac(jg)
end do
!
! Expand the rank 2-->3 and pass values to XBeach variables
call handlearray(sm_ero,s%ero)
call handlearray(sm_depo_ex,s%depo_ex)
call handlearray(sm_Csedeq,s%ceqbg)
!
end subroutine sandmud_erodep
!===================================================================================!
subroutine erosed(morlyr , nmlb , nmub , dt , morfac , &
& ws , umod , h , chezy , taub , &
& nfrac , rhosol , sedd50 , sedd90 , sedtyp , &
& sink , sinkf , sour , sourf , messages , &
& rsedeqb , roller , phase_vel , bed_turb , turb , &
& bore , fric_coeff , wetz , sedcal , urms , &
& par )
!
!! MODIFICATION: rsedeq was added among outputs of erosed
!
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2012.
!
! 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 version 3.
!
! 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.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: erosed.f90 7697 2012-11-16 14:10:17Z boer_aj $
! $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/programs/SandMudBedModule/03_Fortran/example/example/source/erosed.f90 $
!!--description-----------------------------------------------------------------
!
! Function: Computes sedimentation and erosion fluxes
!
!!--declarations----------------------------------------------------------------
use precision
use params
use bedcomposition_module
!use readbedcomposition_module
use message_module
use handlearray_module
use soil_vegetation_module
!
implicit none
!
include 'sedparams.inc'
!
integer , pointer :: flufflyr ! switch for fluff layer concept
real(fp) , dimension(:,:) , pointer :: bfluff0 ! burial coefficient 1 [kg/m2/s]
real(fp) , dimension(:,:) , pointer :: bfluff1 ! burial coefficient 2 [1/s]
real(fp) , dimension(:,:) , pointer :: mfluff ! composition of fluff layer: mass of mud fractions [kg/m2]
!
type(parameters) , intent(in) :: par
type(bedcomp_data) , intent(in) :: morlyr ! bed composition data
integer , intent(in) :: nfrac ! number of sediment fractions
integer , intent(in) :: nmlb ! first cell number
integer , intent(in) :: nmub ! last cell number
integer , dimension(nfrac) , intent(in) :: sedtyp ! sediment type
real(fp) , intent(in) :: dt ! time step [s]
real(fp) , intent(in) :: morfac ! morphological accelaration factor
real(fp) , dimension(nmlb:nmub) , intent(in) :: chezy ! Chezy coefficient for hydraulic roughness [m(1/2)/s]
real(fp) , dimension(nmlb:nmub) , intent(in) :: h ! water depth [m]
real(fp) , dimension(nfrac) , intent(in) :: rhosol ! specific sediment density [kg/m3]
! M: MODIFICATION, sdd50 and sedd50 had dimension(nmlb:nmub)
real(fp) , dimension(nfrac) , intent(in) :: sedd50 ! 50% diameter sediment fraction [m]
real(fp) , dimension(nfrac) , intent(in) :: sedd90 ! 90% diameter sediment fraction [m]
! M: MODIFICATION, sedcal was not present in the orginal version
real(fp) , dimension(nfrac) , intent(in) :: sedcal ! Sediment transport calibration coefficient per grain type [-]
real(fp) , dimension(nmlb:nmub) , intent(in) :: taub ! bottom shear stress [N/m2]
real(fp) , dimension(nmlb:nmub) , intent(in) :: umod ! velocity magnitude (in bottom cell) [m/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(in) :: ws ! sediment settling velocity (hindered) [m/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: sink ! sediment sink flux [m/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: sinkf ! sediment sink flux fluff layer [m/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: sour ! sediment source flux [kg/m2/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: sourf ! sediment source flux fluff layer [kg/m2/s]
type(message_stack) :: messages ! message stack
!
! Local variables
!
character(message_len) :: message ! message
logical :: anymud
integer :: istat
integer :: l ! sediment counter
integer :: nm ! cell counter
real(fp) :: alf1 ! calibration coefficient van Rijn (1984) [-]
real(fp) :: betam ! power factor for adaptation of critical bottom shear stress [-]
real(fp) :: fracf
real(fp) :: mfltot
real(fp) :: rksc ! reference level van Rijn (1984) [m]
real(fp) :: sbot
real(fp) :: smfac ! correction factor for critical bottom shear stress
real(fp) :: ssus
real(fp) , dimension(nmlb:nmub) :: mudcnt
real(fp) , dimension(nmlb:nmub) :: mudfrac ! mud fraction [-]
real(fp) , dimension(nmlb:nmub) :: pmcrit ! critical mud fraction [-]
real(fp) , dimension(nmlb:nmub) :: seddep
real(fp) , dimension(nfrac ) :: E ! erosion velocity [m/s]
real(fp) , dimension(nfrac,nmlb:nmub) :: depeff ! deposition efficiency [-]
real(fp) , dimension(nfrac,nmlb:nmub) :: depfac ! deposition factor (flufflayer=2) [-]
real(fp) , dimension(nfrac,nmlb:nmub) :: eropar ! erosion parameter for mud [kg/m2/s]
real(fp) , dimension(nfrac,nmlb:nmub) :: fixfac ! reduction factor in case of limited sediment availability [-]
! M: changed declaration of frac, before it was 'dimension(nfrac,nmlb:nmub)'
real(fp) , dimension(nmlb:nmub,nfrac) :: frac ! sediment (mass) fraction [-]
real(fp) , dimension(nfrac,nmlb:nmub) :: parfluff0 ! erosion parameter 1 [s/m]
real(fp) , dimension(nfrac,nmlb:nmub) :: parfluff1 ! erosion parameter 2 [ms/kg]
! M: MODIFICATION, rsedeqb was not defined
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: rsedeqb ! equilibrium bed concentration [kg/m3]
real(fp) , dimension(nfrac,nmlb:nmub) :: rsedeq ! equilibrium suspended concentration [kg/m3]
real(fp) , dimension(nfrac,nmlb:nmub) :: tcrdep ! critical bed shear stress for mud sedimentation [N/m2]
real(fp) , dimension(nfrac,nmlb:nmub) :: tcrero ! critical bed shear stress for mud erosion [N/m2]
real(fp) , dimension(nfrac,nmlb:nmub) :: tcrfluff ! critical bed shear stress for fluff layer erosion [N/m2]
! M: MODIFICATION
! Variables employed to relate xbeach original variables to
! sandmud_approach variables in computing sediment transport
! with sb_vr. XBeach variables contain informations on hydrodynamics
real(fp) , dimension(nmlb:nmub) :: roller !
real(fp) , dimension(nmlb:nmub) :: urms !
real(fp) , dimension(nmlb:nmub) :: phase_vel !
real(fp) , dimension(nmlb:nmub) :: bed_turb !
real(fp) , dimension(nmlb:nmub) :: turb !
real(fp) , dimension(nmlb:nmub) :: bore !
real(fp) , dimension(nmlb:nmub) :: fric_coeff !
integer , dimension(nmlb:nmub) :: wetz !
! M: variables to be used for Van Ledden approach with 2 sediment
! fractions, sand and mud
real(fp) :: Mnc ! erosion parameter non-cohesive mixtures
real(fp) :: Mc ! erosion parameter cohesive mixtures
real(fp) :: tau_nc
real(fp) :: tau_c
!
!! executable statements ------------------
!
!
istat = 0
! M: MODIFICATION, at the moment the fluff layer concept is not used
flufflyr => morlyr%settings%flufflyr
bfluff0 => morlyr%state%bfluff0
bfluff1 => morlyr%state%bfluff1
!
! User defined parameters
!
message = 'initializing fluff layer'
if (istat == 0) istat = bedcomp_getpointer_integer(morlyr, 'Flufflyr' , flufflyr)
if (flufflyr>0 .and. istat == 0) istat = bedcomp_getpointer_realfp (morlyr, 'mfluff' , mfluff)
if (flufflyr==1) then
if (istat==0) istat = bedcomp_getpointer_realfp (morlyr, 'Bfluff0' , bfluff0)
if (istat==0) istat = bedcomp_getpointer_realfp (morlyr, 'Bfluff1' , bfluff1)
endif
if (istat /= 0) call adderror(messages, message)
!
call initsed( par , nmlb , nmub , nfrac , flufflyr , &
& alf1 , betam , rksc , pmcrit , bfluff0 , bfluff1 , &
& depeff , depfac , eropar , parfluff0 , parfluff1 , &
& tcrdep , tcrero , tcrfluff )
!
!
mudcnt = 0.0_fp
anymud = .true.
!
! Determine fractions of all sediments in the top layer and compute the mud fraction.
!
call getfrac(morlyr, frac, anymud, mudcnt, mudfrac, nmlb, nmub)
!
! Determine thickness of sediment
!
call getsedthick(morlyr, seddep)
!
! Initialization
!
fixfac = 1.0_fp
rsedeq = 0.0_fp
rsedeqb = 0.0_fp ! M: Added variable to track bed equilibrium concenration
ssus = 0.0_fp
sour = 0.0_fp
sink = 0.0_fp
sinkf = 0.0_fp
sourf = 0.0_fp
!
! Compute change in sediment composition (e.g. based on available fractions and sediment availability)
!
if (nfrac/=2) then
do nm = nmlb, nmub
mfltot = 0.0_fp
if (flufflyr>0) then
do l = 1, nfrac
mfltot = mfltot + mfluff(l,nm)
enddo
endif
do l = 1, nfrac
if (sedtyp(l)==SEDTYP_COHESIVE) then
!
! Compute source and sink fluxes for cohesive sediment (mud)
!
fracf = 0.0_fp
if (mfltot>0.0_fp) fracf = mfluff(l,nm)/mfltot
!
call eromud(ws(l,nm) , fixfac(l,nm) , taub(nm) , frac(nm,l) , fracf , &
& tcrdep(l,nm) , tcrero(l,nm) , eropar(l,nm) , flufflyr , mfltot , &
& tcrfluff(l,nm), depeff(l,nm) , depfac(l,nm) , parfluff0(l,nm), parfluff1(l,nm) , &
& sink(l,nm) , sour(l,nm) , sinkf(l,nm) , sourf(l,nm) )
!
else
!
! Compute correction factor for critical bottom shear stress with sand-mud interaction
!
if ( pmcrit(nm) > 0.0_fp ) then
smfac = ( 1.0_fp + mudfrac(nm) ) ** betam
else
smfac = 1.0_fp
endif
!
! Apply sediment transport formula ( Soulsby-VanrRijn (2007) or vanRijn (1984) )
!
! TODO: here we are inside a loop for the sediment typolgy (ngd or nfrac)
! sb_vr and sednew determine ceqsg and ceqsb for each fraction, then
! it is necessary to create ad hoc subroutine. Temporarily they are used
! in the original version and we take only on of the fractions as output.
!if (true) then
!if (trim(par%form)=='soulsby_vanrijn') then
call sandmud_sb_vr(roller(nm) , phase_vel(nm) , ws(l,nm) , umod(nm) , wetz(nm) , &
& bed_turb(nm) , turb(nm) , bore(nm) , urms(nm) , fric_coeff(nm) , &
& rsedeq(l,nm) , rsedeqb(l,nm) , rhosol(l) , h(nm) , sedd50(l) , &
& sedd90(l) , sedcal(l) , smfac , par )
!elseif (trim(par%form)=='vanthiel_vanrijn') then
! call sandmud_sednew(s,par)
!end if
!
! Trasform suspended sediment concentration from [m3/m3] to [kg/m3]
! to be used in erosand subroutine
!
rsedeq(l,nm) = rsedeq(l,nm)*rhosol(l)
!else
!call vanRijn84(umod(nm) ,sedd50(l) ,sedd90(l) ,h(nm) ,ws(l,nm) , &
! & rhosol(l) ,alf1 ,rksc , &
! & sbot ,ssus ,smfac )
!
!ssus = ssus * rhosol(l)
!
! Compute reference concentration
!
!if (umod(nm)*h(nm)>0.0_fp) then
! rsedeq(l,nm) = frac(l,nm) * ssus / (umod(nm)*h(nm))
!endif
!end if
!
! Compute suspended sediment fluxes for non-cohesive sediment (sand)
!
call erosand(umod(nm) ,chezy(nm) ,ws(l,nm) ,rsedeq(l,nm), &
& sour(l,nm) ,sink(l,nm) ,wetz(nm) )
endif
enddo
enddo
!
! Recompute fluxes due to sand-mud interaction
!
do nm = nmlb, nmub
! Compute erosion velocities
E = 0.0_fp
do l = 1, nfrac
if (frac(nm,l)>0.0_fp) E(l) = sour(l,nm)/(rhosol(l)*frac(nm,l))
enddo
!
! Recompute erosion velocities
!
call sand_mud(nfrac, E, frac(nm,:), mudfrac(nm), sedtyp, pmcrit(nm))
!
! Recompute erosion fluxes
!
do l = 1, nfrac
sour(l,nm) = frac(nm,l)*rhosol(l)*E(l)
enddo
enddo
!
! Add implicit part of source term to sink
!
do l = 1, nfrac
do nm = nmlb, nmub
sink(l,nm) = sink(l,nm)
enddo
enddo
else
!
! Compute sand-mud interaction using original van Ledden approach
!
do l = 1,nfrac
if (sedtyp(l)/=SEDTYP_COHESIVE) then
!
! Critical erosion shear stress for pure sand
!
tcrero(l,:) = par%tausand
end if
end do
!
do nm = nmlb,nmub
!
! Erosion parameter for cohesive and non-cohesive mixtures
!
call erosion_par( mudfrac(nm) , pmcrit(nm) , betam , sedd50 , eropar(:,nm) , sedtyp , &
& Mnc , Mc , tcrero(:,nm) , tau_nc , tau_c , rhosol )
!
if (par%soilveg==1) then
!
! Compute the effect of vegetation in increasing critical shear stress following
! the approach by Van der Meer et al. (2007) and Tuan and Oumeraci (2012)
!
call tau_soil_veg( par, tau_nc, tau_c, frac(nm,:))
!
end if
!
if (mudfrac(nm)>=pmcrit(nm)) then
!
! Choesive regime
!
do l = 1,nfrac
!
if (sedtyp(l)==SEDTYP_COHESIVE) then
!
! Compute source and sink fluxes for cohesive sediment (mud)
!
call C_mud( mudfrac(nm) , pmcrit(nm) , Mc , tau_c , ws(l,nm) , &
& tcrdep(l,nm) , taub(nm) , sour(l,nm) , sink(l,nm) )
!
else
!
call C_sand( ws(l,nm) , mudfrac(nm) , Mc , taub(nm) , &
& tau_c , sour(l,nm) , sink(l,nm) )
!
end if
end do
else
!
! Non-cohesive regime
!
do l = 1,nfrac
!
! Critical shear stress non-cohesive regime
!
if (sedtyp(l)==SEDTYP_COHESIVE) then
!
call NC_mud( mudfrac(nm) , Mnc , taub(nm) , tau_nc , &
& tcrdep(l,nm) , sour(l,nm) , sink(l,nm) , ws(l,nm) )
!
else
!
call NC_sand( ws(l,nm) , mudfrac(nm) , h(nm) , sedd50(l) , rhosol(l) , &
& taub(nm) , tau_nc , sour(l,nm) , sink(l,nm) )
!
end if
end do
end if
end do
end if
!
!
end subroutine erosed
!===================================================================================!
subroutine eromud(ws , fixfac , taub , frac , fracf , &
& tcrdep , tcrero , eropar , flufflyr , mflufftot , &
& tcrfluff , depeff , depfac , parfluff0 , parfluff1 , &
& sink , sour , sinkf , sourf )
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2012.
!
! 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 version 3.
!
! 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.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: eromud.f90 7697 2012-11-16 14:10:17Z boer_aj $
! $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/programs/SandMudBedModule/03_Fortran/example/example/source/eromud.f90 $
!!--description-----------------------------------------------------------------
!
! Function: Computes sediment fluxes at the bed using
! the Partheniades-Krone formulations.
!
!!--declarations----------------------------------------------------------------
use precision
!
implicit none
!
include 'sedparams.inc'
!
integer , intent(in) :: flufflyr ! switch for fluff layer concept
real(fp) , intent(in) :: eropar ! erosion parameter for mud [kg/m2/s]
real(fp) , intent(in) :: depeff ! deposition efficiency [-]
real(fp) , intent(in) :: depfac ! deposition factor (flufflayer=2) [-]
real(fp) , intent(in) :: fixfac ! reduction factor in case of limited sediment availability [-]
real(fp) , intent(in) :: frac ! sediment (mass) fraction [-]
real(fp) , intent(in) :: fracf ! sediment (mass) fraction fluff layer [-]
real(fp) , intent(in) :: mflufftot ! total mass of fluff layer
real(fp) , intent(in) :: parfluff0 ! erosion parameter 1 [s/m]
real(fp) , intent(in) :: parfluff1 ! erosion parameter 2 [s/m]
real(fp) , intent(in) :: taub ! bottom shear stress [N/m2]
real(fp) , intent(in) :: tcrdep ! critical bed shear stress for mud sedimentation [N/m2]
real(fp) , intent(in) :: tcrero ! critical bed shear stress for mud erosion [N/m2]
real(fp) , intent(in) :: tcrfluff ! critical bed shear stress for fluff layer erosion [N/m2]
real(fp) , intent(in) :: ws ! settling velocity [m/s]
real(fp) , intent(out) :: sink ! sediment sink flux [m/s]
real(fp) , intent(out) :: sinkf ! sediment sink flux [m/s]
real(fp) , intent(out) :: sour ! sediment source flux [kg/m2/s]
real(fp) , intent(out) :: sourf ! sediment source flux fluff layer [kg/m2/s]
!
! Local variables
!
real(fp) :: taum
!
!! executable statements ------------------
!
sour = 0.0_fp
sourf = 0.0_fp
sink = 0.0_fp
sinkf = 0.0_fp
!
! Default Partheniades-Krone formula
!
taum = 0.0_fp
if (tcrero>0.0_fp) then
taum = max(0.0_fp, taub/tcrero - 1.0_fp)
endif
sour = eropar * taum
if (tcrdep > 0.0) then
sink = max(0.0_fp , 1.0_fp-taub/tcrdep)
endif
!
! Erosion and deposition to fluff layer
!
if (flufflyr>0) then
taum = max(0.0_fp, taub - tcrfluff)
sourf = min(mflufftot*parfluff1,parfluff0)*taum
sinkf = depeff
sink = 0.0_fp
endif
if (flufflyr==2) then
sinkf = (1.0_fp - depfac)*sinkf
sink = depeff*depfac
endif
!
! Sediment source and sink fluxes
!
sink = ws * sink
sinkf = ws * sinkf
sour = fixfac * frac * sour
sourf = fracf * sourf
!
end subroutine eromud
!===================================================================================!
subroutine vanRijn84(utot ,d50 ,d90 ,h ,ws , &
& rhosol ,alf1 ,rksc , &
& sbot ,ssus ,smfac )
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2011.
!
! 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 version 3.
!
! 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.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: vanRijn84.f90 7697 2012-11-16 14:10:17Z boer_aj $
! $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/programs/SandMudBedModule/03_Fortran/example/example/source/vanRijn84.f90 $
!!--description-----------------------------------------------------------------
!
! computes sediment transport according to
! van rijn (1984)
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
implicit none
!
! Global variables
!
real(fp) , intent(in) :: alf1
real(fp) , intent(in) :: d50 ! grain size diameter (first specified diameter)
real(fp) , intent(in) :: d90 ! grain size diameter (first specified diameter)
real(fp) , intent(in) :: h ! water depth
real(fp) , intent(in) :: rhosol ! density of sediment
real(fp) , intent(in) :: rksc
real(fp) , intent(out) :: sbot ! bed load transport
real(fp) , intent(in) :: smfac ! factor for sand-mud interaction
real(fp) , intent(out) :: ssus ! suspended sediment transport
real(fp) , intent(in) :: utot ! flow velocity
real(fp) , intent(in) :: ws ! settling velocity
!
! Local variables
!
real(fp) :: a
real(fp) :: ah
real(fp) :: beta ! lowest level of integration interval over vertical
real(fp) :: ca
real(fp) :: del
real(fp) :: dstar
real(fp) :: fc
real(fp) :: ff ! coriolis coefficient
real(fp) :: ag ! gravity acceleration
real(fp) :: psi
real(fp) :: rhowat ! density of water
real(fp) :: rmuc
real(fp) :: rnu ! laminar viscosity of water
real(fp) :: t ! time in seconds
real(fp) :: tbc
real(fp) :: tbce
real(fp) :: tbcr
real(fp) :: thetcr
real(fp) :: ustar
real(fp) :: vonkar
real(fp) :: zc
!! MODIFICATION: variable shld (function output) is not declared
!! otherwise lead to error LNK2019 unresolved external symbol _SHLD in ...
!real(fp) :: shld
!
!! executable statements -------------------------------------------------------
!
sbot = 0.0
ssus = 0.0
!
ag = 9.81_fp
rhowat = 1000.0_fp
del = (rhosol - rhowat) / rhowat
rnu = 1e-6_fp
vonkar = 0.41_fp
!
if (h/rksc<1.33 .or. utot<1.E-3) then
return
endif
!
a = rksc
dstar = d50*(del*ag/rnu/rnu)**(1./3.)
!
rmuc = (log10(12.*h/rksc)/log10(12.*h/3./d90))**2
fc = .24*(log10(12.*h/rksc))**( - 2)
tbc = .125*rhowat*fc*utot**2
tbce = rmuc*tbc
thetcr = shld(dstar)
tbcr = (rhosol - rhowat)*ag*d50*thetcr*smfac
t = (tbce - tbcr)/tbcr
!
if (t<.000001) t = .000001
ca = .015*alf1*d50/a*t**1.5/dstar**.3
!
ustar = sqrt(.125*fc)*utot
zc = 0.
beta = 1. + 2.*(ws/ustar)**2
beta = min(beta, 1.5_fp)
psi = 2.5*(ws/ustar)**0.8*(ca/0.65)**0.4
if (ustar>0.) zc = ws/vonkar/ustar/beta + psi
if (zc>20.) zc = 20.
ah = a/h
fc = 0.
if (abs(zc - 1.2)>1.E-4) then
fc = (ah**zc - ah**1.2)/(1. - ah)**zc/(1.2 - zc)
else
fc = -(ah/(1. - ah))**1.2*log(ah)
endif
ff = fc
ssus = ff*utot*h*ca
!
if (t<3.) then
sbot = 0.053*(del)**0.5*sqrt(ag)*d50**1.5*dstar**( - 0.3)*t**2.1
else
sbot = 0.100*(del)**0.5*sqrt(ag)*d50**1.5*dstar**( - 0.3)*t**1.5
endif
end subroutine vanRijn84
!===================================================================================!
subroutine erosand(umod , chezy , ws , rsedeq , &
& sour , sink , wetz )
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2012.
!
! 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 version 3.
!
! 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.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: erosand.f90 7697 2012-11-16 14:10:17Z boer_aj $
! $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/programs/SandMudBedModule/03_Fortran/example/example/source/erosand.f90 $
!!--description-----------------------------------------------------------------
!
! Function: Computes the sour and sink terms for the 2D case
! (Gallappatti aproach)
!
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
!
implicit none
!
! The following list of pointer parameters is used to point inside the gdp structure
!
!
! Global variables
!
real(fp), intent(in) :: chezy ! Chezy coefficient for hydraulic roughness [m(1/2)/s]
real(fp), intent(in) :: umod ! velocity magnitude (in bottom cell) [m/s]
real(fp), intent(in) :: ws ! sediment settling velocity (hindered) [m/s]
real(fp), intent(in) :: rsedeq ! equilibrium concentration [kg/m3]
! M: modification, it is used wetz to control fluxes
integer, intent(in) :: wetz ! Mask for wet(1)/dry(0) cell [-]
real(fp), intent(out) :: sour ! sediment source flux [m/s]
real(fp), intent(out) :: sink ! sediment sink flux [m/s]
!
! Local variables
!
real(fp) :: b
real(fp) :: eps
real(fp) :: hots
real(fp) :: sg
real(fp) :: tsd
real(fp) :: u
real(fp) :: ulog
real(fp) :: ustarc
real(fp) :: w
real(fp) :: wsl
real(fp) :: x
real(fp) :: x2
real(fp) :: x3
!
!! executable statements -------------------------------------------------------
!
eps = 1.0e-6
sg = sqrt(9.81_fp)
!
sour = 0.0_fp
sink = 0.0_fp
!
! local bed shear stress due to currents
!
ustarc = umod*sg/chezy
!
wsl = max(1.0e-3_fp,ws)
if (umod > eps .and. ustarc > eps) then
!
! compute relaxation time using the Gallappatti formulations
!
u = ustarc/umod
!
! limit u to prevent overflow in tsd below
!
u = min(u, 0.15_fp)
if (ustarc > wsl) then
w = wsl/ustarc
else
w = 1.0
endif
b = 1.0
x = w/b
x2 = x*x
x3 = x2*x
ulog = log(u)
tsd = x*exp(( 1.547 - 20.12*u )*x3 &
& + (326.832 *u**2.2047 - 0.2 )*x2 &
& + ( 0.1385*ulog - 6.4061)*x &
& + ( 0.5467*u + 2.1963) )
!
hots = wsl/tsd
sour = rsedeq*hots*wetz
sink = hots*wetz
else
sink = wsl*wetz
endif
end subroutine erosand
!===================================================================================!
subroutine sand_mud(nfrac, E, frac, mudfrac, sedtyp, pmcrit)
!
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2012.
!
! 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 version 3.
!
! 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.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: sand_mud.f90 7697 2012-11-16 14:10:17Z boer_aj $
! $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/programs/SandMudBedModule/03_Fortran/example/example/source/sand_mud.f90 $
!!--description-----------------------------------------------------------------
!
! Function: Computes erosion velocities based
! on sand-mud interaction (Van Ledden (2003), Van Kessel (2002))
! Array E is recomputed.
! Method used:
!
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
!
implicit none
!
include 'sedparams.inc'
!
integer , intent(in) :: nfrac ! number of sediment fractions
integer , dimension(nfrac) , intent(in) :: sedtyp ! sediment type
real(fp) , dimension(nfrac) , intent(in) :: frac ! sediment (mass) fraction [-]
real(fp) , intent(in) :: mudfrac ! mud fraction [-]
real(fp) , intent(in) :: pmcrit ! critical mud fraction [-]
real(fp) , dimension(nfrac) , intent(inout) :: E ! sediment erosion velocity [m/s]
!
! Local variables
!
integer :: istat ! error flag
integer :: l ! sediment counter
real(fp) :: Es_avg ! average erosion velocity for sand fractions [m/s]
real(fp) :: Em_avg ! average erosion velocity for mud fractions [m/s]
!
!
!! executable statements ------------------
!
!
! No sand mud interaction if there is no mud, only mud or pmcrit<0
if (pmcrit<0.0_fp) return
if (comparereal(mudfrac,0.0_fp)==0) return
if (comparereal(mudfrac,1.0_fp)==0) return
!
Es_avg = 0.0_fp
Em_avg = 0.0_fp
!
! Compute average erosion velocity for sand fractions
!
do l = 1, nfrac
if (sedtyp(l)/= SEDTYP_COHESIVE) then
Es_avg = Es_avg + frac(l)*E(l)
endif
enddo
Es_avg = Es_avg/(1-mudfrac)
!
if ( mudfrac <= pmcrit ) then
!
! Non-cohesive regime
! (mud is proportionally eroded with the sand)
!
do l = 1, nfrac
if (sedtyp(l) == SEDTYP_COHESIVE ) then
if (Es_avg>0.0_fp) then
E(l) = Es_avg
else
E(l) = 0.0_fp
endif
endif
enddo
else
!
! Cohesive regime
!
! erosion velocity for mud is interpolated between the non-cohesive and fully mud regime
! fully mud regime : mudfrac = 1 -> E(l) is not changed
! non-cohesive regime: mudfrac = pmcrit -> E(l) = Es_avg
!
do l = 1, nfrac
if ( sedtyp(l)==SEDTYP_COHESIVE ) then
if (Es_avg>0.0_fp .and. E(l)>0.0_fp) then
E(l) = E(l)*(Es_avg/E(l))**((1.0_fp-mudfrac)/(1.0_fp-pmcrit))
else
E(l) = 0.0_fp
endif
Em_avg = Em_avg + frac(l)*E(l)
endif
enddo
Em_avg = Em_avg/mudfrac
!
! sand is proportionally eroded with the mud
!
do l = 1, nfrac
if (sedtyp(l) /= SEDTYP_COHESIVE) then
E(l) = Em_avg
endif
enddo
endif
end subroutine sand_mud
!===================================================================================!
subroutine initsed( par , nmlb , nmub , nfrac , flufflyr , &
& alf1 , betam , rksc , pmcrit , bfluff0 , bfluff1 , &
& depeff , depfac , eropar , parfluff0 , parfluff1 , &
& tcrdep , tcrero , tcrfluff )
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2012.
!
! 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 version 3.
!
! 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.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: initsed.f90 7697 2012-11-16 14:10:17Z boer_aj $
! $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/programs/SandMudBedModule/03_Fortran/example/example/source/initsed.f90 $
!!--description-----------------------------------------------------------------
!
! Function: Settings of morphological parameters
!
!!--declarations----------------------------------------------------------------
use precision
use params
!
implicit none
!
!
! Global variables
!
type(parameters) , intent(in) :: par
integer , intent(in) :: flufflyr ! switch for fluff layer concept
integer , intent(in) :: nfrac ! number of sediment fractions
integer , intent(in) :: nmlb ! first cell number
integer , intent(in) :: nmub ! first cell number
real(fp) , intent(out) :: alf1 ! calibration coefficient van Rijn (1984) [-]
real(fp) , intent(out) :: betam ! power factor for adaptation of critical bottom shear stress [-]
real(fp) , intent(out) :: rksc ! reference level van Rijn (1984) [m]
real(fp) , dimension(nmlb:nmub) , intent(out) :: pmcrit ! critical mud fraction [-]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: bfluff0 ! burial coefficient 1 [kg/m2/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: bfluff1 ! burial coefficient 2 [1/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: depeff ! deposition efficiency [-]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: depfac ! deposition factor (flufflayer=2) [-]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: eropar ! erosion parameter for mud [kg/m2/s]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: parfluff0 ! erosion parameter 1 [s/m]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: parfluff1 ! erosion parameter 2 [ms/kg]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: tcrdep ! critical bed shear stress for mud sedimentation [N/m2]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: tcrero ! critical bed shear stress for mud erosion [N/m2]
real(fp) , dimension(nfrac,nmlb:nmub) , intent(out) :: tcrfluff ! critical bed shear stress for fluff layer erosion [N/m2]
!
!! executable statements ------------------
!
! ================================================================================
! USER INPUT
! ================================================================================
!
! Parameters sediment
!
! M: modification to take erosion parameters from params.txt
eropar = par%Mmud ! erosion parameter for mud [kg/m2/s]
tcrdep = par%taudep ! critical bed shear stress for mud sedimentation [N/m2]
tcrero = par%taumud ! critical bed shear stress for mud erosion [N/m2]
!
! Parameters fluff layer
!
depeff = 0.95_fp ! deposition efficiency [-]
depfac = 0.2_fp ! deposition factor (flufflayer=2) [-]
parfluff0 = 2.0e-1_fp ! erosion parameter 1 [s/m]
parfluff1 = 1.0_fp ! erosion parameter 2 [ms/kg]
tcrfluff = 0.05_fp ! critical bed shear stress for fluff layer erosion [N/m2]
if (flufflyr==1) then
bfluff0 = 0.0_fp ! burial coefficient 1 [kg/m2/s]
bfluff1 = 0.0_fp ! burial coefficient 2 [1/s]
endif
!
! Parameters sand-mud interaction
!
betam = 1.0_fp ! power factor for adaptation of critical bottom shear stress [-]
pmcrit = par%pmcr ! critical mud fraction [-]
!
! Parameters sediment transport formulation
!
alf1 = 2.0_fp ! calibration coefficient [-]
rksc = 0.1_fp ! reference level [m]
!
! ================================================================================
end subroutine initsed
!===================================================================================!
function shld(dstar)
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2011.
!
! 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 version 3.
!
! 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.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
!
! contact: delft3d.support@deltares.nl
! Stichting Deltares
! P.O. Box 177
! 2600 MH Delft, The Netherlands
!
! All indications and logos of, and references to, "Delft3D" and "Deltares"
! are registered trademarks of Stichting Deltares, and remain the property of
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: shld.f90 7697 2012-11-16 14:10:17Z boer_aj $
! $HeadURL: https://svn.oss.deltares.nl/repos/openearthtools/trunk/programs/SandMudBedModule/03_Fortran/example/example/source/shld.f90 $
!!--description-----------------------------------------------------------------
!
! determines shields parameter according
! to shields curve
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
implicit none
!
! Global variables
!
real(fp), intent(in) :: dstar ! critical dimensionless grain size parameter
real(fp) :: shld
!
!! executable statements -------------------------------------------------------
!
if (dstar<=4.) then
shld = 0.240/dstar
elseif (dstar<=10.) then
shld = 0.140/dstar**0.64
elseif (dstar<=20.) then
shld = 0.040/dstar**0.10
elseif (dstar<=150.) then
shld = 0.013*dstar**0.29
else
shld = 0.055
endif
end function shld
!===================================================================================!
subroutine settling(s,par,ws,rhos)
use params
use spaceparams
implicit none
type(parameters) :: par
type(spacepars) :: s
real*8, dimension(par%ngd) , intent(in) :: rhos
real*8, dimension(par%ngd,par%nx+1,par%ny+1), intent(out) :: ws
real*8 :: sd
real*8 :: nu
integer :: jg
!
nu = 1.0040E-6
do jg = 1,par%ngd
sd = (rhos(jg)-par%rho)/par%rho
if (par%d50(jg)<65.0E-6) then
ws(jg,:,:) = 2.5E-4
elseif (par%d50(jg)>65.0E-6 .and. par%d50(jg)<100.0E-6) then
ws(jg,:,:) = (sd-1)*par%g*par%d50(jg)**2/(18.d0*nu)
elseif (par%d50(jg)>=100.0E-6 .and. par%d50(jg)<1000.0E-6) then
ws(jg,:,:) = 10.d0*nu/par%d50(jg)*((1+0.01d0*(sd-1)*par%g*par%d50(jg)**3.d0/nu**2.d0)**0.5d0-1)
else
ws(jg,:,:) = 1.1d0*((sd-1)*par%g*par%d50(jg))**0.5d0
end if
end do
!
end subroutine settling
!===================================================================================!
subroutine tau_bed(s,par,tau_tot,sm_hh,sm_U,sm_urms)
use params
use spaceparams
use precision
use handlearray_module
use bedcomposition_module
!
implicit none
type(parameters) :: par
type(spacepars) :: s
type(bedcomp_data) :: bc
real*8, dimension((par%nx+1)*(par%ny+1)), intent(out) :: tau_tot ! [N/m2] total bed shear stress
real*8, dimension((par%nx+1)*(par%ny+1)), intent(in) :: sm_hh ! [m] water depth
real*8, dimension((par%nx+1)*(par%ny+1)), intent(in) :: sm_U ! [m/s] eulerian velocity magnitude
real*8, dimension((par%nx+1)*(par%ny+1)) :: sm_urms ! [m/s] root mean squared orbital velocity
! local
real*8, dimension((par%nx+1)*(par%ny+1)) :: fw ! [-] wave friction parameter
real*8, dimension((par%nx+1)*(par%ny+1)) :: A ! [m] semi-orbital excursion
real*8, dimension((par%nx+1)*(par%ny+1)) :: Cd ! [-] drag coefficient
real*8, dimension((par%nx+1)*(par%ny+1)) :: tau_w ! [N/m2] wave bed shear stress
real*8, dimension((par%nx+1)*(par%ny+1)) :: tau_c ! [N/m2] current bed shear stress
real*8, dimension((par%nx+1)*(par%ny+1)) :: tau_m ! [N/m2] mean shear stress due to wave-current interaction
real*8 :: phi ! [rad] angle between wave direction and current direction
real*8 :: eq_d50 ! [m] equivalent grain roughness to compute the friction factor
real*8, dimension((par%nx+1)*(par%ny+1)) :: hloc
integer :: cellnumber,pos
!
!
hloc = max(sm_hh,0.01)
phi = 0.0_fp
eq_d50 = 5.25E-6
cellnumber = (par%nx+1)*(par%ny+1)
!
!A = sm_urms*par%Trep/(2.0_fp*par%px)
!fw = 1.39_fp*(A/par%z0)**(-0.52_fp)
!
do pos=1,cellnumber
if (sm_urms(pos)<=0.00001_fp) then
fw(pos) = 0.0_fp
else
fw(pos) = 1.39_fp*(sm_urms(pos)*par%Trep/(2*par%px*eq_d50))**(-0.52_fp)
end if
end do
!
tau_w = 0.5_fp*par%rho*fw*sm_urms**2
Cd = (0.4_fp/(log(max(hloc,10.d0*par%z0)/par%z0)-1.0_fp))**2
tau_c = par%rho*Cd*sm_U**2
!
do pos=1,cellnumber
if (tau_w(pos)+tau_c(pos)==0.0_fp) then
tau_m(pos) = 0.0_fp
else
tau_m(pos) = tau_c(pos)*(1.0_fp+1.2_fp*(tau_w(pos)/(tau_w(pos)+tau_c(pos)))**3.2_fp)
end if
end do
!
tau_tot = sqrt((tau_m+tau_w*cos(phi))**2+(tau_w*sin(phi))**2)
end subroutine tau_bed
!===================================================================================!
subroutine sandmud_sb_vr(roller , phase_vel , ws , umod , wetz , &
& bed_turb , turb , bore , urms , fric_coeff , &
& eqcs , eqcb , rhosoil , depth , myd50 , &
& myd90 , sedcal , smfac , par )
!use params
implicit none
type(parameters) :: par
real*8 , intent(in) :: roller
real*8 , intent(in) :: phase_vel
real*8 , intent(in) :: ws
real*8 , intent(in) :: umod
real*8 :: bed_turb
real*8 , intent(in) :: turb
real*8 , intent(in) :: bore
real*8 , intent(in) :: urms
real*8 , intent(in) :: fric_coeff
real*8 , intent(in) :: rhosoil
real*8 , intent(in) :: depth
real*8 , intent(in) :: myd50
real*8 , intent(in) :: myd90
real*8 , intent(in) :: sedcal
integer, intent(in) :: wetz
real*8 , intent(in) :: smfac
real*8 , intent(out) :: eqcs
real*8 , intent(out) :: eqcb
! local
real*8 :: ML
real*8 :: hloc
real*8 :: dcfin
real*8 :: dcf
real*8 :: vmg
real*8 :: urms2
real*8 :: Ts
real*8 :: Ucr
real*8 :: Cd
real*8 :: Asb
real*8 :: Ass
real*8 :: delta
real*8 :: term1
real*8 :: term2
real*8 :: dster
!
hloc = max(depth,0.01)
vmg = 0.d0
delta = (rhosoil-par%rho)/par%rho
dster = (delta*par%g/1.d-12)**(1.d0/3.d0)*myd50
!
! Wave breaking induced turbulence due to short waves
! Compute mixing length
!
ML = sqrt(2*roller*par%Trep/(rhosoil*phase_vel))
ML = min(ML,hloc)
dcfin = exp(min(100.d0,hloc/max(ML,0.01d0)))
dcf = min(1.d0,1.d0/(dcfin-1.d0))
!
! Turbulence
!
bed_turb = turb*dcf
if (trim(par%turb)=='bore_averaged') then
bed_turb = bed_turb*par%Trep/bore
end if
!
! Including long wave stirring
!
if (par%lws==1) then
vmg = umod
elseif (par%lws==0) then
vmg = (1.d0-1.d0/par%cats/par%Trep*par%dt)*vmg + (1.d0/par%cats/par%Trep*par%dt)*umod
end if
!
urms2 = urms**2+1.45d0*bed_turb
!
Ts = par%tsfac*hloc/ws
Ts = max(Ts,par%Tsmin)
!
! Determination of critical velocity
!
if (myd50<=0.0005d0) then
Ucr = 0.19d0*myd50**0.1d0*log10(4.d0*hloc/myd90)
else if (myd50<=0.05d0) then
Ucr = 8.5d0*myd50**0.6d0*log10(4.d0*hloc/myd90)
else
Ucr = 10
write(*,'(a,i4)') 'In process'
end if
!
! Account for the presence of mud in affecting critical shear stress
! in case of non-cohesive sand-mud mixtures
!
Ucr = Ucr*sqrt(smfac)
!
! Drag coefficient
!
Cd = (0.40d0/(log(max(hloc,10.d0*par%z0)/par%z0)-1.0d0))**2
!
! Transport parameters
!
Asb = 0.005d0*hloc*(myd50/hloc/(delta*par%g*myd50))**1.2d0 ! bed load coefficent
Ass = 0.012d0*myd50*dster**(-0.6d0)/(delta*par%g*myd50)**1.2d0 ! suspended load coeffient
term1 = (vmg**2+0.018/Cd*par%sws*urms2)
term1 = min(term1,par%smax*par%g/fric_coeff*myd50*delta)
term1 = sqrt(term1)
term2 = 0.d0
if (term1>Ucr .and. depth>par%eps) then
term2 = (term1-Ucr)**2.4d0
end if
!
! Deterination of equilibrium sediment concentrations
!
eqcb = Asb*term2
eqcb = min(eqcb/hloc,par%cmax/2) ! maximum equilibrium bed concentration
eqcb = (1-par%bulk)*eqcb*sedcal*wetz
eqcs = Ass*term2
eqcs = min(eqcs/hloc,par%cmax/2) ! maximum equilibrium suspended concentration
eqcs = (eqcs+par%bulk*eqcs)*sedcal*wetz
!
end subroutine sandmud_sb_vr
!===================================================================================!
subroutine erosion_par( pm , pmcr , betam , sedd50 , eropar , sedtyp , &
& Mnc , Mc , tcr , tau_nc , tau_c , rhosol )
use precision
implicit none
include 'sedparams.inc'
real*8, intent(in) :: pm ! [-] actual mud fraction
real*8, intent(in) :: pmcr ! [-] critical mud fraction
real*8, dimension(2), intent(in) :: sedd50 ! [m] sediment diameter
real*8, dimension(2), intent(in) :: tcr ! [N/m2] critical shear stress
real*8, dimension(2), intent(in) :: eropar ! [] erosion parameter
real*8, dimension(2), intent(in) :: rhosol ! [] soil density
integer, dimension(2), intent(in) :: sedtyp ! sediment type
real*8, intent(out) :: Mnc ! [] erosion parameter non-cohesive regime
real*8, intent(out) :: Mc ! [] erosion parameter cohesive regime
real*8, intent(out) :: tau_nc ! [N/m2] critical shear stress non-cohesive mixtures
real*8, intent(out) :: tau_c ! [N/m2] critical shear stress cohesive mixtures
real*8 :: betam
real*8 :: alfa
real*8 :: nu ! [m2/s] kinematic viscosity
real*8 :: myd50
real*8 :: mytcr_s
real*8 :: mytcr_m
real*8 :: rho_sand
real*8 :: rho_mud
real*8 :: Me ! [kg/m2/s] erosion parameter pure mud
real*8 :: Dstar
real*8 :: expon
real*8 :: rhow
real*8 :: sr
real*8 :: g
integer :: l
!
!
do l = 1,2
if (sedtyp(l)/=SEDTYP_COHESIVE) then
myd50 = sedd50(l)
mytcr_s = tcr(l)
rho_sand = rhosol(l)
else
mytcr_m = tcr(l)
Me = eropar(l)
rho_mud = rhosol(l)
end if
end do
!
g = 9.806d0
rhow = 1025.d0
sr = rho_sand/rhow
alfa = 0.075d0
nu = 1.0040E-6
Dstar = myd50*((sr-1)*g/nu**2)**(1.d0/3.d0)
Mnc = alfa*sqrt((sr-1)*g*myd50)/(Dstar**0.9)
Mnc = 0.05d0
!
expon = (1-pm)/(1-pmcr)
Mc = Me*(Mnc/Me/(1-pmcr))**expon
!
!
tau_c = (mytcr_s*(1+pmcr)**betam-mytcr_m)/(1-pmcr)*(1-pm)+mytcr_m
tau_nc = mytcr_s*(1+pm)**betam
!
end subroutine erosion_par
!===================================================================================!
subroutine C_mud( pm , pmcr , Mc , tau_cr , ws , &
& tcrdep , taub , sour , sink )
implicit none
real*8, intent(in) :: pm ! [-] actual mud fraction
real*8, intent(in) :: pmcr ! [-] critical mud fraction
real*8, intent(in) :: Mc ! [] erosion parameter for cohesive mixtures
real*8, intent(in) :: tau_cr ! [N/m2] critical shear stress
real*8, intent(in) :: tcrdep ! [N/m2] critical shear stress for deposition
real*8, intent(in) :: taub ! [N/m2] bottom shear stress
real*8, intent(in) :: ws ! [N/m2] bottom shear stress
real*8, intent(out) :: sour ! [kg/m2/s] erosion flux
real*8, intent(out) :: sink ! [m/s] deposition rate
real*8 :: taum
!
!
sour = 0.d0
sink = 0.d0
!
taum = 0.d0
if (tau_cr>0.d0) then
taum = max(0.d0,taub/tau_cr-1)
end if
!
! Erosion flux for mud cohesive regime
!
sour = pm*Mc*taum
!
if (tcrdep>0.d0) then
sink = max(0.d0,1-taub/tcrdep)
end if
!
sink = ws*sink
!
end subroutine C_mud
!===================================================================================!
subroutine C_sand( ws , pm , Mc , taub , &
& tau_cr , sour , sink )
implicit none
real*8, intent(in) :: ws ! [m/s] sand settling veocity
real*8, intent(in) :: pm ! [-] actual mud fraction
real*8, intent(in) :: Mc ! [] erosion parameter cohesive mixtures
real*8, intent(in) :: taub ! [N/m2] bed shear stress
real*8, intent(in) :: tau_cr ! [N/m2] critical shear stress
real*8, intent(out) :: sour ! [kg/m2/s] erosion flux
real*8, intent(out) :: sink ! [m/s] deposition rate
real*8 :: taum
!
!
taum = 0.d0
if (tau_cr>0.d0) then
taum = max(0.d0,taub/tau_cr-1)
end if
!
! Sand is proportionally eroded with mud in cohesive regime
!
sour = (1-pm)*Mc*taum
!
! TODO: fix this, choosing a formulation for sand settling
sink = ws*1
!
end subroutine C_sand
!===================================================================================!
subroutine NC_mud( pm , Mnc , taub , tau_cr , &
& tcrdep , sour , sink , ws )
!
! FIX THIS: SHOULDN'T MUD BE ERODED PROPORTIONALLY TO SAND IN
! NON COHESIVE REGIME ?
!
implicit none
real*8, intent(in) :: pm
real*8, intent(in) :: Mnc
real*8, intent(in) :: taub
real*8, intent(in) :: tau_cr
real*8, intent(in) :: tcrdep
real*8, intent(in) :: ws
real*8, intent(out) :: sour
real*8, intent(out) :: sink
real*8 :: taum
!
!
taum = 0.d0
if (tau_cr>0.d0) then
taum = max(0.d0,taub/tau_cr-1)
end if
!
! Erosion flux for mud non-cohesive regime
!
sour = pm/(1-pm)*Mnc*taum
!
if (tcrdep>0.d0) then
sink = max(0.d0,1-taub/tcrdep)
end if
!
sink = sink*ws
!
end subroutine NC_mud
!===================================================================================!
subroutine NC_sand( ws , pm , h , myd50 , rhosand , &
& taub , tau_cr , sour , sink )
implicit none
real*8, intent(in) :: ws
real*8, intent(in) :: pm
real*8, intent(in) :: h
real*8, intent(in) :: myd50
real*8, intent(in) :: rhosand
real*8, intent(in) :: taub
real*8, intent(in) :: tau_cr
real*8, intent(out) :: sour
real*8, intent(out) :: sink
real*8 :: taum
real*8 :: Dstar
real*8 :: alfa
real*8 :: nu
real*8 :: g
real*8 :: sr
!
!
g = 9.806
sr = rhosand/g
alfa = 1.0E-5
nu = 1.0040E-6
Dstar = myd50*((sr-1)*g/nu**2)*(1.d0/3.d0)
!
taum = 0.d0
if (tau_cr>0.d0) then
taum = max(0.d0,taub/tau_cr-1)
end if
!
! Erosion flux for sand non-cohesive regime
!
sour = (1-pm)*ws*1.5d0*(myd50/h/Dstar**0.3)*taum**1.5d0
!
! TODO: fix this, choosing a formulation for sand settling
sink = ws*1
!
end subroutine NC_sand
!===================================================================================!
end module sandmud_module