!----- AGPL --------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2015.
!
! This file is part of Delft3D (D-Flow Flexible Mesh component).
!
! Delft3D is free software: you can redistribute it and/or modify
! it under the terms of the GNU Affero General Public License as
! published by the Free Software Foundation version 3.
!
! Delft3D 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 Affero General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License
! along with Delft3D. 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",
! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting
! Deltares, and remain the property of Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
! $Id: fm_erosed.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fm_erosed.f90 $
subroutine fm_erosed()
!! (nmmax ,kmax ,icx ,icy ,lundia , &
!! & nst ,lsed ,lsedtot ,lsal ,ltem , &
!! & lsecfl ,kfs ,kfu ,kfv ,dzs1 , &
!! & r0 ,u0eul ,v0eul ,s0 ,dps , &
!! & z0urou ,z0vrou ,sour ,sink ,rhowat , &
!! & ws ,rsedeq ,z0ucur ,z0vcur ,sigmol , &
!! & taubmx ,s1 ,uorb ,tp ,sigdif , &
!! & lstsci ,thick ,dicww ,kcs , &
!! & kcu ,kcv ,guv ,gvu ,sbuu , &
!! & sbvv ,seddif ,hrms ,ltur , &
!! & teta ,rlabda ,aks ,saleqs , &
!! & sbuut ,sbvvt ,entr ,wstau ,hu , &
!! & hv ,rca ,dss ,ubot ,rtur0 , &
!! & temeqs ,gsqs ,guu ,gvv ,kfsmin , &
!! & kfsmax ,dzs0 ,kfumin ,kfumax ,kfvmin , &
!! & kfvmax ,dzu1 ,dzv1 ,gdp )
!!!!--description-----------------------------------------------------------------
!!!
!!! Function: Computes sediment fluxes at the bed using
!!! the Partheniades-Krone formulations.
!!! Arrays SOURSE and SINKSE are filled and added
!!! to arrays SOUR and SINK
!!! Computes bed load transport for sand sediment
!!! Arrays SBUU and SBVV are filled.
!!! Computes vertical sediment diffusion coefficient
!!! Array SEDDIF is filled
!!! Includes wave asymmetry effects on sand bed-load
!!! transport
!!! Bed slope effects computed at the U and V velocity
!!! points
!!! Method used: Attention: pointer ll for 'standard' FLOW
!!! arrays is shifted with lstart
!!!
!!!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
use mathconsts, only: pi
use bedcomposition_module
use morphology_data_module
use sediment_basics_module
use m_physcoef, only: ag, vonkar, sag, ee
use m_sediment, only: stmpar, sedtra, stm_included, mtd
use m_flowgeom, only: ndxi, bl, kfs, lnxi, lnx, ln, dxi, Ndx !, csu, snu
use m_flow, only: s0, s1, ucx, ucy, kbot, ktop, kmx, kmxn, plotlin
use m_flowtimes, only: julrefdat, dts, time1
use unstruc_files, only: mdia
use message_module, only: write_error
use m_transport, only: ised1
use dfparall
use m_alloc
use m_missing
use m_physcoef, only: frcuni, ifrctypuni
use m_turbulence, only: vicwws
!
implicit none
!
logical, pointer :: wave
real(fp) :: eps = 1.0e-6_fp
! real(fp) :: ag from m_physcoef
real(fp) :: vicmol = 1.3e-6_fp
real(fp) :: gammax = 1.0_fp
logical :: scour = .false.
real(fp) , dimension(:) , pointer :: rksr
! real(fp) :: vonkar from m_physcoef
!! real(fp) , pointer :: timsec
!! real(fp) , pointer :: timhr
!! integer , pointer :: julday
logical :: ubot_from_com = .true. !! promoted approach, so only option in FM
!! logical , pointer :: flmd2l
logical :: flmd2l = .false.
! sedpar
integer , pointer :: nmudfrac
real(fp) , dimension(:) , pointer :: rhosol
real(fp) , dimension(:) , pointer :: cdryb
real(fp) , dimension(:,:,:) , pointer :: logseddia
real(fp) , dimension(:) , pointer :: logsedsig
real(fp) , dimension(:) , pointer :: sedd10
real(fp) , dimension(:) , pointer :: sedd50
real(fp) , dimension(:) , pointer :: sedd90
real(fp) , dimension(:) , pointer :: sedd50fld
real(fp) , dimension(:) , pointer :: dstar
real(fp) , dimension(:) , pointer :: taucr
real(fp) , dimension(:) , pointer :: tetacr
real(fp) , dimension(:) , pointer :: mudcnt
real(fp) , dimension(:) , pointer :: pmcrit
integer , dimension(:) , pointer :: nseddia
integer , dimension(:) , pointer :: sedtyp
logical , pointer :: anymud
real(fp) , dimension(:) , pointer :: sedtrcfac
logical , pointer :: bsskin
real(fp) , dimension(:) , pointer :: thcmud
real(fp) , pointer :: kssilt
real(fp) , pointer :: kssand
! morpar
real(fp) , pointer :: thresh
real(fp) , pointer :: sus
real(fp) , pointer :: bed
real(fp) , pointer :: susw
real(fp) , pointer :: sedthr
real(fp) , pointer :: bedw
integer , pointer :: i10
integer , pointer :: i50
integer , pointer :: i90
integer , pointer :: nxx
real(fp) , dimension(:) , pointer :: xx
real(fp) , pointer :: morfac
logical , pointer :: varyingmorfac
logical , pointer :: multi
real(fp) , pointer :: factcr
integer , pointer :: ihidexp
real(fp) , pointer :: asklhe
real(fp) , pointer :: mwwjhe
real(fp) , pointer :: ffthresh
real(fp) , pointer :: espir
logical , pointer :: epspar
real(fp) , pointer :: camax
real(fp) , pointer :: aksfac
real(fp) , pointer :: rdc
integer , pointer :: iopkcw
logical , pointer :: oldmudfrac
integer , pointer :: iflufflyr
real(fp) , dimension(:,:) , pointer :: sinkf
real(fp) , dimension(:,:) , pointer :: sourf
real(fp) , dimension(:,:) , pointer :: depfac
real(fp) , dimension(:,:) , pointer :: mfluff
! trapar
integer , dimension(:) , pointer :: iform
real(fp) , dimension(:,:) , pointer :: par
integer , pointer :: max_integers
integer , pointer :: max_reals
integer , pointer :: max_strings
character(256) , dimension(:) , pointer :: dll_function
integer(pntrsize), dimension(:) , pointer :: dll_handle
integer , dimension(:) , pointer :: dll_integers
real(hp) , dimension(:) , pointer :: dll_reals
character(256) , dimension(:) , pointer :: dll_strings
character(256) , dimension(:) , pointer :: dll_usrfil
! sedtra
real(fp) , dimension(:) , pointer :: bc_mor_array
real(fp) , dimension(:,:) , pointer :: dbodsd
real(fp) , dimension(:) , pointer :: dcwwlc
real(fp) , dimension(:) , pointer :: dm
real(fp) , dimension(:) , pointer :: dg
real(fp) , dimension(:) , pointer :: dgsd
real(fp) , dimension(:,:) , pointer :: dxx
real(fp) , dimension(:) , pointer :: e_dzdn
real(fp) , dimension(:) , pointer :: e_dzdt
real(fp) , dimension(:) , pointer :: epsclc
real(fp) , dimension(:) , pointer :: epswlc
real(fp) , dimension(:,:) , pointer :: fixfac
real(fp) , dimension(:,:) , pointer :: frac
integer , dimension(:) , pointer :: kfsed
integer , dimension(:,:) , pointer :: kmxsed
real(fp) , dimension(:) , pointer :: mudfrac
real(fp) , dimension(:) , pointer :: sandfrac
real(fp) , dimension(:,:) , pointer :: hidexp
real(fp) , dimension(:) , pointer :: rsdqlc
real(fp) , dimension(:,:) , pointer :: sbcx
real(fp) , dimension(:,:) , pointer :: sbcy
real(fp) , dimension(:,:) , pointer :: e_sbcn
real(fp) , dimension(:,:) , pointer :: e_sbct
real(fp) , dimension(:,:) , pointer :: sbwx
real(fp) , dimension(:,:) , pointer :: sbwy
real(fp) , dimension(:,:) , pointer :: e_sbwn
real(fp) , dimension(:,:) , pointer :: e_sbwt
real(fp) , dimension(:) , pointer :: sddflc
real(fp) , dimension(:,:) , pointer :: srcmax
real(fp) , dimension(:,:) , pointer :: sswx
real(fp) , dimension(:,:) , pointer :: sswy
real(fp) , dimension(:,:) , pointer :: e_sswn
real(fp) , dimension(:,:) , pointer :: e_sswt
real(fp) , dimension(:,:) , pointer :: sxtot
real(fp) , dimension(:,:) , pointer :: sytot
real(fp) , dimension(:,:) , pointer :: sinkse
real(fp) , dimension(:,:) , pointer :: sourse
real(fp) , dimension(:,:) , pointer :: sour_im
real(fp) , dimension(:,:) , pointer :: taurat
real(fp) , dimension(:) , pointer :: ust2
real(fp) , dimension(:) , pointer :: umod
real(fp) , dimension(:) , pointer :: uuu
real(fp) , dimension(:) , pointer :: vvv
real(fp) , dimension(:) , pointer :: wslc
real(fp) , dimension(:) , pointer :: zumod
!
integer , pointer :: iunderlyr
real(prec) , dimension(:,:) , pointer :: bodsed
!! include 'flow_steps_f.inc'
!
! Local parameters
!
integer, parameter :: kmax2d = 20
!
! Global variables
!
!! integer , intent(in) :: icx ! Increment in the X-dir., if ICX= NMAX then computation proceeds in the X-dir. If icx=1 then computation proceeds in the Y-dir.
!! integer , intent(in) :: icy ! Increment in the Y-dir. (see ICX)
!! integer , intent(in) :: kmax ! Description and declaration in esm_alloc_int.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfsmin ! Description and declaration in iidim.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfsmax ! Description and declaration in iidim.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfumin ! Description and declaration in iidim.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfumax ! Description and declaration in iidim.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfvmin ! Description and declaration in iidim.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfvmax ! Description and declaration in iidim.f90
!! integer , intent(in) :: lsal ! Description and declaration in dimens.igs
integer, pointer :: lsed
integer, pointer :: lsedtot
!! integer , intent(in) :: lstsci ! Description and declaration in esm_alloc_int.f90
!! integer , intent(in) :: ltem ! Description and declaration in dimens.igs
integer :: lsecfl = 0
integer :: ltur = 0
!! integer :: lundia ! Description and declaration in inout.igs
!! integer , intent(in) :: nmmax ! Description and declaration in dimens.igs
!! integer , intent(in) :: nst
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kcs ! Description and declaration in esm_alloc_int.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kcu ! Description and declaration in esm_alloc_int.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kcv ! Description and declaration in esm_alloc_int.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfs ! Description and declaration in esm_alloc_int.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfu ! Description and declaration in esm_alloc_int.f90
!! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfv ! Description and declaration in esm_alloc_int.f90
!! real(fp) , intent(in) :: dt
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, lsed) :: aks ! Description and declaration in esm_alloc_real.f90
!! real(prec), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: dps ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: entr ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: guv ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gvu ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gsqs ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: guu ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gvv ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(:), pointer :: hrms
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: hu ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: hv ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(:), pointer :: rlabda
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax, ltur), intent(in) :: rtur0 ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: s0 ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: s1 ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: sbuut
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: sbvvt
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: taubmx ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(:), pointer :: teta
real(fp), dimension(:), pointer :: tp
real(fp), dimension(:), pointer :: uorb
real(fp), dimension(:), pointer :: ubot
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: wstau ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: z0ucur
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: z0urou ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: z0vcur
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: z0vrou ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax) , intent(in) :: dicww ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax, *) , intent(in) :: ws ! Description and declaration in esm_alloc_real.f90
real(fp) , dimension(:,:), pointer :: seddif
real(fp) , dimension(:), pointer :: rhowat
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: u0eul ! EULARIAN U-velocities at old time level
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: v0eul ! EULARIAN V-velocities at old time level
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, *) :: r0 ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, *) :: sink ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, *) :: sour ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lsed) :: rsedeq ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, lsed) :: dss ! Description and declaration in esm_alloc_real.f90
real(fp) , dimension(:,:), pointer :: caksrho
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, lsedtot) :: sbuu ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, lsedtot)
!! real(fp) , dimension(kmax) , intent(in) :: thick ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(lstsci) , intent(out) :: sigdif ! Description and declaration in esm_alloc_real.f90
!! real(fp) , dimension(lstsci) :: sigmol ! Description and declaration in esm_alloc_real.f90
!! real(fp) , intent(in) :: saleqs
!! real(fp) , intent(in) :: temeqs
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: dzs1 ! Description and declaration in rjdim.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: dzs0 ! Description and declaration in rjdim.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: dzu1 ! Description and declaration in rjdim.f90
!! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: dzv1 ! Description and declaration in rjdim.f90
!
! Local variables
!
integer :: i
integer :: iln
integer :: istat
integer :: j
integer :: k
integer :: k2d
integer :: kbed
integer :: kmaxsd
!! integer :: kn
!! integer :: ku
!! integer :: kv
integer :: l
!! integer :: ll
integer :: lstart
!! integer :: m
!! integer :: n
!! integer :: ndm
!! integer :: nhystp
integer :: nm
!! integer :: nmd
!! integer :: nm_pos ! indicating the array to be exchanged has nm index at the 2nd place, e.g., dbodsd(lsedtot,nm)
!! integer :: nmu
!! integer :: num
logical :: error
integer :: klc
integer :: kmaxlc
logical :: suspfrac ! suspended component sedtyp(l)/=SEDTYP_NONCOHESIVE_TOTALLOAD
real(fp) :: aks_ss3d
real(fp) :: caks
real(fp) :: caks_ss3d
real(fp) :: chezy
!! real(fp) :: conc2d
real(fp) :: delr
real(fp) :: di50
!! real(fp) :: difbot
real(fp) :: drho
!! real(fp) :: dstari
real(fp) :: dtmor
real(fp) :: dzdx
real(fp) :: dzdy
!! real(fp) :: fi
real(fp) :: fracf
real(fp) :: grkg
real(fp) :: grm2
real(fp) :: grlyrs
real(fp) :: h0
real(fp) :: h1
real(fp) :: rc
real(fp) :: mfltot
real(fp) :: salinity
!! real(fp) :: sinkfluff
real(fp) :: sinktot
real(fp) :: sourfluff
real(fp) :: spirint ! local variable for spiral flow intensity r0(nm,1,lsecfl)
real(fp) :: taks
real(fp) :: taks0
real(fp) :: tauadd
real(fp) :: taub
real(fp) :: tauc
real(fp) :: tdss ! temporary variable for dss
real(fp) :: temperature
real(fp), dimension(kmx+1) :: thicklc
real(fp), dimension(kmx+1) :: siglc
real(fp) :: thick0
real(fp) :: thick1
real(fp) :: trsedeq ! temporary variable for rsedeq
real(fp) :: tsd
real(fp) :: tsigmol ! temporary variable for sigmol
real(fp) :: twsk
real(fp) :: u
real(fp) :: ubed
real(fp) :: umean
real(fp) :: ustarc
real(fp) :: utot
real(fp) :: v
real(fp) :: vbed
real(fp) :: velb
real(fp) :: velm
real(fp) :: vmean
real(fp) :: z0cur
real(fp) :: z0rou
real(fp) :: zvelb
real(fp), dimension(:), allocatable :: E ! erosion velocity [m/s]
real(fp), dimension(0:kmax2d) :: dcww2d
real(fp), dimension(0:kmax2d) :: sddf2d
real(fp), dimension(0:kmax2d) :: ws2d
real(fp), dimension(kmax2d) :: rsdq2d
real(fp), dimension(kmax2d), save :: sig2d = &
(/ -0.0874, -0.2472, -0.3797, -0.4897, -0.5809, -0.6565, -0.7193, &
& -0.7713, -0.8145, -0.8503, -0.8800, -0.9046, -0.9250, -0.9419, -0.9560,&
& -0.9676, -0.9773, -0.9854, -0.9920, -0.9975 /)
real(fp), dimension(kmax2d), save :: thck2d = &
(/ 0.1747, 0.1449, 0.1202, 0.0997, 0.0827, 0.0686, 0.0569, 0.0472, &
& 0.0391, 0.0325, 0.0269, 0.0223, 0.0185, 0.0154, 0.0127, 0.0106, 0.0088,&
& 0.0073, 0.0060, 0.0050 /)
!! real(fp), dimension(kmax) :: concin3d
real(fp), dimension(kmax2d) :: concin2d
character(256) :: errmsg
double precision :: cz_dum
integer :: idummy = 0
!
!
!! executable statements -------------------------------------------------------
!
! exit the routine immediately if sediment transport (and morphology) is NOT included in the simulation
!
error = .false.
if (.not.stm_included) return
wave => mtd%have_waves
!! vicmol => gdp%gdphysco%vicmol
!! gammax => gdp%gdnumeco%gammax
!! eps => gdp%gdconst%eps
!! scour => gdp%gdscour%scour
hrms => mtd%hrms
tp => mtd%tp
teta => mtd%teta
rlabda => mtd%rlabda
uorb => mtd%uorb
ubot => mtd%ubot
rksr => mtd%rksr
rhowat => mtd%rhowat
seddif => mtd%seddif
caksrho => mtd%caksrho
!! timsec => gdp%gdinttim%timsec
!! timhr => gdp%gdinttim%timhr
!! julday => gdp%gdinttim%julday
!! ubot_from_com => gdp%gdprocs%ubot_from_com
!! flmd2l => gdp%gdprocs%flmd2l
! stmpar
lsed => stmpar%lsedsus
lsedtot => stmpar%lsedtot
! sedpar
nmudfrac => stmpar%sedpar%nmudfrac
rhosol => stmpar%sedpar%rhosol
cdryb => stmpar%sedpar%cdryb
logseddia => stmpar%sedpar%logseddia
logsedsig => stmpar%sedpar%logsedsig
sedd10 => stmpar%sedpar%sedd10
sedd50 => stmpar%sedpar%sedd50
sedd90 => stmpar%sedpar%sedd90
sedd50fld => stmpar%sedpar%sedd50fld
dstar => stmpar%sedpar%dstar
taucr => stmpar%sedpar%taucr
tetacr => stmpar%sedpar%tetacr
mudcnt => stmpar%sedpar%mudcnt
pmcrit => stmpar%sedpar%pmcrit
nseddia => stmpar%sedpar%nseddia
sedtyp => stmpar%sedpar%sedtyp
anymud => stmpar%sedpar%anymud
sedtrcfac => stmpar%sedpar%sedtrcfac
bsskin => stmpar%sedpar%bsskin
thcmud => stmpar%sedpar%thcmud
kssilt => stmpar%sedpar%kssilt
kssand => stmpar%sedpar%kssand
! morpar
thresh => stmpar%morpar%thresh
sus => stmpar%morpar%sus
bed => stmpar%morpar%bed
susw => stmpar%morpar%susw
sedthr => stmpar%morpar%sedthr
bedw => stmpar%morpar%bedw
i10 => stmpar%morpar%i10
i50 => stmpar%morpar%i50
i90 => stmpar%morpar%i90
nxx => stmpar%morpar%nxx
xx => stmpar%morpar%xx
multi => stmpar%morpar%multi
factcr => stmpar%morpar%factcr
ihidexp => stmpar%morpar%ihidexp
asklhe => stmpar%morpar%asklhe
mwwjhe => stmpar%morpar%mwwjhe
ffthresh => stmpar%morpar%thresh
morfac => stmpar%morpar%morfac
varyingmorfac => stmpar%morpar%varyingmorfac
espir => stmpar%morpar%espir
epspar => stmpar%morpar%epspar
camax => stmpar%morpar%camax
aksfac => stmpar%morpar%aksfac
rdc => stmpar%morpar%rdc
iopkcw => stmpar%morpar%iopkcw
oldmudfrac => stmpar%morpar%oldmudfrac
sinkf => stmpar%morpar%flufflyr%sinkf
sourf => stmpar%morpar%flufflyr%sourf
iflufflyr => stmpar%morpar%flufflyr%iflufflyr
depfac => stmpar%morpar%flufflyr%depfac
mfluff => stmpar%morpar%flufflyr%mfluff
! trapar
iform => stmpar%trapar%iform
par => stmpar%trapar%par
max_integers => stmpar%trapar%max_integers
max_reals => stmpar%trapar%max_reals
max_strings => stmpar%trapar%max_strings
dll_function => stmpar%trapar%dll_function
dll_handle => stmpar%trapar%dll_handle
dll_integers => stmpar%trapar%dll_integers
dll_reals => stmpar%trapar%dll_reals
dll_strings => stmpar%trapar%dll_strings
dll_usrfil => stmpar%trapar%dll_usrfil
! sedtra
bc_mor_array => sedtra%bc_mor_array
dbodsd => sedtra%dbodsd
dcwwlc => sedtra%dcwwlc
dm => sedtra%dm
dg => sedtra%dg
dgsd => sedtra%dgsd
dxx => sedtra%dxx
e_dzdn => sedtra%e_dzdn
e_dzdt => sedtra%e_dzdt
epsclc => sedtra%epsclc
epswlc => sedtra%epswlc
fixfac => sedtra%fixfac
frac => sedtra%frac
kfsed => sedtra%kfsed
kmxsed => sedtra%kmxsed
mudfrac => sedtra%mudfrac
sandfrac => sedtra%sandfrac
hidexp => sedtra%hidexp
rsdqlc => sedtra%rsdqlc
sbcx => sedtra%sbcx
sbcy => sedtra%sbcy
e_sbcn => sedtra%e_sbcn
e_sbct => sedtra%e_sbct
sbwx => sedtra%sbwx
sbwy => sedtra%sbwy
e_sbwn => sedtra%e_sbwn
e_sbwt => sedtra%e_sbwt
sddflc => sedtra%sddflc
sswx => sedtra%sswx
sswy => sedtra%sswy
e_sswn => sedtra%e_sswn
e_sswt => sedtra%e_sswt
sxtot => sedtra%sxtot
sytot => sedtra%sytot
sinkse => sedtra%sinkse
sourse => sedtra%sourse
sour_im => sedtra%sour_im
srcmax => sedtra%srcmax
taurat => sedtra%taurat
ust2 => sedtra%ust2
umod => sedtra%umod
uuu => sedtra%uuu
vvv => sedtra%vvv
wslc => sedtra%wslc
zumod => sedtra%zumod
!
if (varyingmorfac) then
call updmorfac(stmpar%morpar, time1/3600.0_fp, julrefdat)
endif
!! !
!! nm_pos = 1
if (scour) then
!
! Second parameter is zero: save taubmx(*) in gdp%gdscour
!
!! call shearx(taubmx, 0, gdp)
endif
!
! Determine total thickness of the mud layers
! to be used in computation of skin friction (Soulsby 2004)
!
if (bsskin) then
call detthcmud(stmpar%morlyr, thcmud)
endif
!
! Initialisation:
! reset sediment sources and sinks
! set default kmxsed layer
! set kfsed
!
lstart = ised1 - 1 ! = max(isalt, itemp)
!
! Reset Sourse and Sinkse arrays for all (l,nm)
!
sinkse = 0.0_fp
sourse = 0.0_fp
sour_im = 0.0_fp
! source and sink terms fluff layer
if (iflufflyr>0) then
sinkf = 0.0_fp
sourf = 0.0_fp
endif
!
! Reset Sediment diffusion arrays for (l,nmk)
!
seddif = 0.0_fp
caksrho = 0.0_fp
!
! Reset Bed Shear Ratio for all nm and l = 1:lsedtot
!
taurat = 0.0_fp
!
! Set zero bedload transport for all nm and l = 1:lsedtot
!
!! sbuu = 0.0_fp
!! sbvv = 0.
sbcx = 0.0_fp
sbcy = 0.0_fp
e_sbcn = 0.0_fp
e_sbct = 0.0_fp
sbwx = 0.0_fp
sbwy = 0.0_fp
e_sbwn = 0.0_fp
e_sbwt = 0.0_fp
sswx = 0.0_fp
sswy = 0.0_fp
e_sswn = 0.0_fp
e_sswt = 0.0_fp
sxtot = 0.0_fp
sytot = 0.0_fp
!! !
!! call dfexchg( dps,1, 1, dfloat, nm_pos, gdp)
!! !
do nm = 1, ndxi
if ((s1(nm) - bl(nm))*kfs(nm) > sedthr) then
kfsed(nm) = 1
else
kfsed(nm) = 0
endif
enddo
!! !
!! call dfexchg( kfsed,1, 1, dfint, nm_pos, gdp)
!
! Determine fractions of all sediments the top layer and
! compute the mud fraction.
!
if (lsedtot > 1) then
call getfrac(stmpar%morlyr,frac ,anymud ,mudcnt , &
& mudfrac ,1 ,ndxi)
endif
! BEGIN DEBUG
frac = 1d0
! END DEBUG
!! !
!! ! Calculate velocity components and magnitude at the zeta points
!! ! based on velocity in the bottom computational layer
!! !
!! call z_dwnvel(nmmax ,kmax ,icx ,kcs ,kfu , &
!! & kfv ,kcu ,kcv ,s1 ,dps , &
!! & u0eul ,v0eul ,uuu ,vvv ,umod , &
!! & zumod ,dzs1 ,hu ,hv ,kfsed , &
!! & kfsmin ,kfsmax ,kfumin ,kfumax ,kfvmin , &
!! & kfvmax ,dzu1 ,dzv1 ,z0ucur ,z0vcur , &
!! & vonkar ,gdp )
! Quick hack for 2D ...
! compute ucx, ucy
call setucxucyucxuucyu()
uuu = ucx
vvv = ucy
umod = sqrt(uuu*uuu + vvv*vvv)
zumod = (s1 - bl)/ee
!! call dfexchg( uuu, 1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( vvv, 1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( umod, 1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( zumod,1, 1, dfloat, nm_pos, gdp)
!
! Get the reduction factor if thickness of sediment at bed is less than
! user specified threshold. Also get maximum erosion source SRCMAX
! (used for cohesive sediments).
!
dtmor = dts * morfac
!
call getfixfac(stmpar%morlyr, 1 , ndxi , lsedtot, &
& ndxi , fixfac , ffthresh )
!
! Set fixfac to 1.0 for tracer sediments and adjust frac
!
istat = bedcomp_getpointer_integer(stmpar%morlyr, 'IUnderLyr', iunderlyr)
if (ffthresh>0.0_hp .or. iunderlyr/=1) then
srcmax = 1.0e+10_fp
elseif (iunderlyr==1) then
istat = bedcomp_getpointer_realprec(stmpar%morlyr,'bodsed',bodsed)
do l = 1, lsed
if (ffthresh<1.0e-10_fp) then
!
! Compute SRCMAX (only used for cohesive sediments)
!
do nm = 1, ndxi
!
! If user-specified THRESH is <= 0.0, the erosion flux is effectively not limited by FIXFAC since ffthresh is 1e-10
! but by the amount of sediment that is available
!
srcmax(nm, l) = bodsed(l, nm)*cdryb(l)/dtmor
enddo
endif
!
if (sedtrcfac(l)>0.0_fp) then
grkg = 1.0_fp / (rhosol(l)*pi*sedd50(l)**3/6.0_fp) ! Number of grains per kg
grm2 = 0.5_fp / (pi*sedd50(l)**2) ! Number of grains per m^2 -- Not quite correct: maximum area of grain is pi*r^2 not pi*d^2, using porosity factor of 0.5
do nm = 1, ndxi
fixfac(nm, l) = 1.0_fp
grlyrs = bodsed(l, nm) * grkg / grm2 ! Number of grain layers
frac(nm, l) = min(max(0.0_fp, grlyrs), 1.0_fp)*sedtrcfac(l)
enddo
endif
enddo
endif
!
! in case of multiple (non-mud) fractions, the following quantities
! --- that are initialized in INISED --- may be time-dependent and
! they must be updated here or after updating the bed levels in
! BOTT3D. Since we do it here, these quantities will lag a half time
! step behind on the output files. If these statements are moved to
! BOTT3D, the GETFRAC call above must be shifted too.
!
if (lsedtot-nmudfrac > 1) then
!
! calculate arithmetic mean sediment diameter Dm
! calculate geometric mean sediment diameter Dg
! calculate percentiles Dxx
!
call compdiam(frac ,sedd50 ,sedd50 ,sedtyp ,lsedtot , &
& logsedsig ,nseddia ,logseddia ,ndxi ,1 , &
& ndxi ,xx ,nxx ,sedd50fld ,dm , &
& dg ,dxx ,dgsd )
!
! determine hiding & exposure factors
!
call comphidexp(frac ,dm ,ndxi ,lsedtot , &
& sedd50 ,hidexp ,ihidexp ,asklhe , &
& mwwjhe ,1 ,ndxi )
!
! compute sand fraction
!
call compsandfrac(frac ,sedd50 ,ndxi ,lsedtot , &
& sedtyp ,sandfrac ,sedd50fld , &
& 1 ,ndxi )
endif
!
do iln = 1, lnxi
!
! compute normal component of bed slopes at edges
! SPvdP: Ask Bert about orientation of db/dn
!
e_dzdn(iln) = (bl(ln(1,iln)) - bl(ln(2,iln)))*dxi(iln)
enddo
!! !
!! call dfexchg( dzduu,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( dzdvv,1, 1, dfloat, nm_pos, gdp)
!
! Start of main loop over sediment fractions for suspended sediment
! sources, sinks, equilibrium concentrations and vertical diffusion
! coefficients, and bed-load transport vector components at water
! level points
!
!! call dfexchg( z0ucur,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( z0vcur,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( z0urou,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( z0vrou,1, 1, dfloat, nm_pos, gdp)
!! do l = 1, lsedtot
!! call dfexchg( ws(:,:,l),0, kmax, dfloat, nm_pos, gdp)
!! enddo
!
do nm = 1, ndxi
if (kfs(nm)/=1) cycle
!
! do not calculate sediment sources, sinks, and bed load
! transport in areas with very shallow water.
!
if (kfsed(nm) == 0) then
!
! Very shallow water:
! set sediment diffusion coefficient
! and set zero equilibrium concentrations
!
if (kmx>1) then
! at layer interfaces, but not at bed and surface
do k = kbot(nm), ktop(nm)-1
do l = 1, lsed
seddif(l, k) = 0.7_fp * vicwws(k) ! sigdifi * vicwws(k) + difsed(l)
enddo
enddo
! in layers
do k = kbot(nm), ktop(nm)
do l = 1, lsed
!! rsedeq(l, k) = 0.0_fp
enddo
enddo
endif
cycle
endif
!
! kfsed(nm) == 1
!
h0 = max(0.01_fp, s0(nm) - bl(nm))
h1 = max(0.01_fp, s1(nm) - bl(nm))
! nmd = nm - icx
!! ndm = nm - icy
!! call nm_to_n_and_m(nm, n, m, gdp)
! In z-layer use kmaxlc, siglc and thicklc for morphology (like sigma-layer).
kmaxlc = kmx
if (kmx>0) then
!! !
!! ! 3D CASE
!! !
!! kbed = kbot(nm)
!! thicklc = 0.0_fp
!! klc = 1
!! do k = kfsmax(nm),kfsmin(nm),-1
!! thicklc(klc) = dzs1(nm, k)/h1
!! klc=klc+1
!! enddo
!! siglc = 0.0_fp
!! kmaxlc =klc-1
!! siglc(1) = -0.5_fp*thicklc(1)
!! do klc = 2, kmaxlc
!! siglc(klc) = siglc(klc - 1) - 0.5_fp*(thicklc(klc) + thicklc(klc - 1))
!! enddo
else
kbed = nm
kmaxlc = 1
thicklc = 1.0_fp
endif
! BEGIN DEBUG
kbed = nm
! END DEBUG
!
! Compute depth-averaged velocity components at cell centre
!
umean = ucx(nm)
vmean = ucy(nm)
velm = sqrt(umean**2+vmean**2)
!
ubed = ucx(kbed)
vbed = ucy(kbed)
velb = sqrt(ubed**2 + vbed**2)
if (kmaxlc>1) then
zvelb = 0.5_fp*thicklc(kmaxlc)*h1
else
zvelb = h1/ee
endif
!
! Calculate current related roughness
!
!! kn = max(1, kfu(nm) + kfu(nmd) + kfv(nm) + kfv(ndm))
!! z0cur = ( kfu(nmd)*z0ucur(nmd) + kfu(nm)*z0ucur(nm) &
!! & + kfv(ndm)*z0vcur(ndm) + kfv(nm)*z0vcur(nm) )/kn
z0cur = 0.01_fp
! SPvdP: use z0 from unstruc
call getczz0 (h1, frcuni, ifrctypuni, cz_dum, z0cur)
!
! Calculate total (possibly wave enhanced) roughness
!
!! z0rou = ( kfu(nmd)*z0urou(nmd) + kfu(nm)*z0urou(nm) &
!! & + kfv(ndm)*z0vrou(ndm) + kfv(nm)*z0vrou(nm) )/kn
z0rou = z0cur
chezy = sag * log( 1.0_fp + h1/max(1.0e-8_fp,ee*z0rou) ) / vonkar
!
! bed shear stress as used in flow, or
! skin fiction following Soulsby; "Bed shear stress under
! combined waves and currents on rough and smoooth beds"
! Estproc report TR137, 2004
!
if (bsskin) then
!
! Compute bed stress resulting from skin friction
!
call compbsskin (umean , vmean , h1 , wave , &
& uorb(nm), tp (nm) , teta(nm), kssilt , &
& kssand , thcmud(nm), taub , rhowat(kbed), &
& vicmol )
else
!
! use max bed shear stress, rather than mean
!
!! taub = taubmx(nm)
ustarc = umod(nm)*vonkar/log(1.0_fp + zumod(nm)/z0rou)
taub = ustarc*ustarc*rhowat(kbed)
endif
!
ustarc = umod(nm)*vonkar/log(1.0_fp + zumod(nm)/z0rou)
if (scour) then
!
! Calculate extra stress (tauadd) for point = nm, if so required by
! user input. Increment TAUB(MX) and USTARC.
!
!! call shearx(tauadd, nm, gdp)
taub = sqrt(taub**2 + tauadd**2)
!
tauc = rhowat(kbed)*ustarc**2
tauc = sqrt(tauc**2 + tauadd**2)
ustarc = sqrt(tauc/rhowat(kbed))
else
tauadd = 0.0_fp
endif
!
! Compute effective depth averaged velocity
!
utot = ustarc * chezy / sag
u = utot * uuu(nm) / (umod(nm)+eps)
v = utot * vvv(nm) / (umod(nm)+eps)
!
!! if (lsal > 0) then
!! salinity = r0(nm, kbed, lsal)
!! else
!! salinity = saleqs
salinity = 0.0_fp
!! endif
!! if (ltem > 0) then
!! temperature = r0(nm, kbed, ltem)
!! else
!! temperature = temeqs
temperature = 15.0_fp
!! endif
!! !
taks0 = 0.0_fp
!
! Calculate Van Rijn's reference height
!
if (iopkcw==1) then
rc = 30.0_fp*z0cur
else
rc = rdc
endif
taks0 = max(aksfac*rc, 0.01_fp*h1)
!
if (wave .and. tp(nm)>0.0_fp) then
delr = 0.025_fp
taks0 = max(0.5_fp*delr, taks0)
endif
!
! Limit maximum aks to 20% of water depth
! (may be used when water depth becomes very small)
!
taks0 = min(taks0, 0.2_fp*h1)
!
! Input parameters are passed via dll_reals/integers/strings-arrays
!
if (max_reals < MAX_RP) then
write(errmsg,'(a)') 'Insufficient space to pass real values to transport routine.'
call write_error(errmsg, unit=mdia)
error = .true.
return
endif
dll_reals(RP_TIME ) = real(time1 ,hp)
dll_reals(RP_EFUMN) = real(u ,hp)
dll_reals(RP_EFVMN) = real(v ,hp)
dll_reals(RP_EFVLM) = real(utot ,hp)
dll_reals(RP_UCHAR) = real(uuu(nm) ,hp)
dll_reals(RP_VCHAR) = real(vvv(nm) ,hp)
dll_reals(RP_VELCH) = real(umod(nm) ,hp)
dll_reals(RP_ZVLCH) = real(zumod(nm) ,hp)
dll_reals(RP_DEPTH) = real(h1 ,hp)
dll_reals(RP_CHEZY) = real(chezy ,hp)
if (wave) then
dll_reals(RP_HRMS ) = real(min(gammax*h1, hrms(nm)) ,hp)
dll_reals(RP_TPEAK) = real(tp(nm) ,hp)
dll_reals(RP_TETA ) = real(teta(nm) ,hp)
dll_reals(RP_RLAMB) = real(rlabda(nm) ,hp)
dll_reals(RP_UORB ) = real(uorb(nm) ,hp)
else
dll_reals(RP_HRMS ) = 0.0_hp
dll_reals(RP_TPEAK) = 0.0_hp
dll_reals(RP_TETA ) = 0.0_hp
dll_reals(RP_RLAMB) = 0.0_hp
dll_reals(RP_UORB ) = 0.0_hp
endif
dll_reals(RP_D10MX) = real(dxx(nm,i10),hp)
dll_reals(RP_D90MX) = real(dxx(nm,i90),hp)
dll_reals(RP_MUDFR) = real(mudfrac(nm),hp)
dll_reals(RP_RHOWT) = real(rhowat(kbed) ,hp) ! Density of water
dll_reals(RP_SALIN) = real(salinity ,hp)
dll_reals(RP_TEMP ) = real(temperature ,hp)
dll_reals(RP_GRAV ) = real(ag ,hp)
dll_reals(RP_VICML) = real(vicmol ,hp)
dll_reals(RP_TAUB ) = real(taub ,hp) !taubmx incremented with tauadd
dll_reals(RP_UBED ) = real(ubed ,hp)
dll_reals(RP_VBED ) = real(vbed ,hp)
dll_reals(RP_VELBD) = real(velb ,hp)
dll_reals(RP_ZVLBD) = real(zvelb ,hp)
dll_reals(RP_VNKAR) = real(vonkar ,hp)
dll_reals(RP_Z0CUR) = real(z0cur ,hp)
dll_reals(RP_Z0ROU) = real(z0rou ,hp)
dll_reals(RP_DG ) = real(dg(nm) ,hp)
dll_reals(RP_SNDFR) = real(sandfrac(nm) ,hp)
dll_reals(RP_DGSD ) = real(dgsd(nm) ,hp)
if (ltur >= 1) then
!! dll_reals(RP_KTUR ) = real(rtur0(nm,kbed,1),hp)
endif
dll_reals(RP_UMEAN) = real(umean ,hp)
dll_reals(RP_VMEAN) = real(vmean ,hp)
dll_reals(RP_VELMN) = real(velm ,hp)
dll_reals(RP_USTAR) = real(ustarc ,hp)
!
if (max_integers < MAX_IP) then
write(errmsg,'(a)') 'Insufficient space to pass integer values to transport routine.'
call write_error(errmsg, unit=mdia)
error = .true.
return
endif
dll_integers(IP_NM ) = nm
!! dll_integers(IP_N ) = n
!! dll_integers(IP_M ) = m
!
if (max_strings < MAX_SP) then
write(errmsg,'(a)') 'Insufficient space to pass strings to transport routine.'
call write_error(errmsg, unit=mdia)
error = .true.
return
endif
!! dll_strings(SP_RUNID) = gdp%runid
!
! total mass in fluff layer
!
mfltot = 0.0_fp
if (iflufflyr>0) then
do l = 1, lsedtot
mfltot = mfltot + max(0.0_fp,mfluff(l,nm))
enddo
endif
!
do l = 1, lsedtot
!
! fraction specific quantities
!
dll_reals(RP_HIDEX) = real(hidexp(nm,l) ,hp)
dll_reals(RP_RHOSL) = real(rhosol(l) ,hp)
dll_integers(IP_ISED ) = l
dll_strings(SP_USRFL) = dll_usrfil(l)
!
do i = 1,stmpar%trapar%npar
j = stmpar%trapar%iparfld(i,l)
if (j>0) then
par(i,l) = stmpar%trapar%parfld(nm,j)
endif
enddo
!
if (sedtyp(l) == SEDTYP_COHESIVE) then
!
! sediment type COHESIVE
!
dll_reals(RP_D50 ) = 0.0_hp
dll_reals(RP_DSS ) = 0.0_hp
dll_reals(RP_DSTAR) = 0.0_hp
!! dll_reals(RP_SETVL) = real(ws(nm, kbed, l) ,hp) ! Vertical velocity near bedlevel
!! if (flmd2l) then
!! par(11,l) = entr(nm)
!! endif
!
klc = 0
dcwwlc = 0.0_fp
wslc = 0.0_fp
!! do k = kfsmax(nm),kfsmin(nm)-1,-1
!! dcwwlc(klc) = dicww(nm, k)
!! wslc(klc) = ws(nm, k, l)
!! klc=klc+1
!! enddo
!
! Fluff layer parameters
!
fracf = 0.0_fp
if (iflufflyr>0) then
if (mfltot>0.0_fp) fracf = max(0.0_fp,mfluff(l,nm))/mfltot
endif
!
kmaxsd = kmaxlc ! for mud fractions kmaxsd points to the grid cell at the bottom of the water column
thick0 = thicklc(kmaxsd) * h0
thick1 = thicklc(kmaxsd) * h1
! following lines to be replaced by call ...
sinktot = 0.0_fp
sourfluff = 0.0_fp
!! call erosilt(thicklc ,kmaxlc ,wslc ,mdia , &
!! & thick0 ,thick1 ,fixfac(nm,l), srcmax(nm, l),&
!! & frac(nm,l) ,oldmudfrac ,flmd2l ,iform(l) , &
!! & par ,max_integers,max_reals ,max_strings , &
!! & dll_function(l),dll_handle(l),dll_integers,dll_reals, &
!! & dll_strings ,iflufflyr ,mfltot ,fracf , &
!! & error ,wstau(nm) ,sinktot ,sourse(nm,l), sourfluff)
if (error) return
!
if (iflufflyr>0) then
if (iflufflyr==2) then
sinkf(l,nm) = sinktot*(1.0_fp - depfac(l,nm))
sinkse(nm,l) = sinktot*depfac(l,nm)
else
sinkf(l,nm) = sinktot
sinkse(nm,l) = 0.0_fp
endif
!
sourf(l,nm) = sourfluff
else
sinkse(nm,l) = sinktot
sourse(nm,l) = sourse(nm,l) + sourfluff ! sourfluff should actually always be 0 already
endif
!
if (kmx>1) then
!
! For 3D model set sediment diffusion coefficient
! NOTE THAT IF ALGEBRAIC OR K-L TURBULENCE MODEL IS USED THEN WAVES
! ONLY AFFECT THE VERTICAL TURBULENT MIXING VIA THE ENHANCED BED
! ROUGHNESS
!
klc = 0
!! do k = kfsmax(nm),kfsmin(nm)-1,-1
!! seddif(nm, k, l) = dcwwlc(klc)
!! klc=klc+1
!! enddo
endif
!
! l runs from 1 to lsedtot, kmxsed is defined for 1:lsed
! The first lsed fractions are the suspended fractions (including cohesive ones),
! so this goes right
!
!! kmxsed(nm, l) = kfsmin(nm)+kmaxlc-kmaxsd
cycle
endif
!
! sediment type NONCOHESIVE_SUSPENDED or NONCOHESIVE_TOTALLOAD
!
!! ll = lstart + l
suspfrac = sedtyp(l)/=SEDTYP_NONCOHESIVE_TOTALLOAD
!
! Calculation for sand or bedload
!
! Reset Prandtl-Schmidt number for sand fractions
!
if (suspfrac) then
!! sigdif(ll) = 1.0_fp
endif
tsd = -999.0_fp
di50 = sedd50(l)
!DEBUG if (di50 < 0.0_fp) then
!
! Space varying sedd50 specified in array sedd50fld:
! Recalculate dstar, tetacr and taucr for each nm,l - point
! This code is copied from inised (uniform sedd50)
!
!DEBUG di50 = sedd50fld(nm)
drho = (rhosol(l)-rhowat(kbed)) / rhowat(kbed)
dstar(l) = di50 * (drho*ag/vicmol**2)**0.3333_fp
if (dstar(l) < 1.0_fp) then
if (iform(l) == -2) then
tetacr(l) = 0.115_fp / (dstar(l)**0.5_fp)
else
tetacr(l) = 0.24_fp / dstar(l)
endif
elseif (dstar(l) <= 4.0_fp) then
if (iform(l) == -2) then
tetacr(l) = 0.115_fp / (dstar(l)**0.5_fp)
else
tetacr(l) = 0.24_fp / dstar(l)
endif
elseif (dstar(l)>4.0_fp .and. dstar(l)<=10.0_fp) then
tetacr(l) = 0.14_fp / (dstar(l)**0.64_fp)
elseif (dstar(l)>10.0_fp .and. dstar(l)<=20.0_fp) then
tetacr(l) = 0.04_fp / (dstar(l)**0.1_fp)
elseif (dstar(l)>20.0_fp .and. dstar(l)<=150.0_fp) then
tetacr(l) = 0.013_fp * (dstar(l)**0.29_fp)
else
tetacr(l) = 0.055_fp
endif
taucr(l) = factcr * (rhosol(l)-rhowat(kbed)) * ag * di50 * tetacr(l)
!DEBUG endif
!
if (suspfrac) then
!! tsigmol = sigmol(ll)
!! tdss = dss(nm, l)
!! twsk = ws(nm, kbed, l)
! BEGIN DEBUG
tsigmol = 1.0_fp
tdss = di50
twsk = 1d-4
! END DEBUG
else
!
! use dummy values for bedload fractions
!
tsigmol = 1.0_fp
tdss = di50
twsk = 0.0_fp
endif
!
! NONCOHESIVE fraction specific quantities
!
dll_reals(RP_D50 ) = real(di50 ,hp)
dll_reals(RP_DSS ) = real(tdss ,hp)
dll_reals(RP_DSTAR) = real(dstar(l),hp)
dll_reals(RP_SETVL) = real(twsk ,hp) ! Vertical velocity near bedlevel
par(1,l) = ag
par(2,l) = rhowat(kbed) ! rhow
par(3,l) = rhosol(l)
par(4,l) = (rhosol(l)-rhowat(kbed)) / rhowat(kbed)
par(5,l) = 1.0E-6 ! rnu from md-tran.*
par(6,l) = di50
!
! TODO: compute bed slope vector in cell centre from normal components of bed slope at edges
!
dzdx = 0.0_fp
dzdy = 0.0_fp
!
! SWITCH 2DH/3D SIMULATIONS
!
if (kmaxlc > 1) then
!! !
!! ! 3D CASE
!! !
!! if (suspfrac) then
!! !
!! ! Fill local 1dv arrays with fall velocity and diffusivity.
!! !
!! klc = 0
!! dcwwlc = 0.0_fp
!! wslc = 0.0_fp
!! do k = kfsmax(nm),kfsmin(nm)-1,-1
!! dcwwlc(klc) = dicww(nm, k)
!! wslc(klc) = ws(nm, k, l)
!! klc=klc+1
!! enddo
!! endif
!! taks = 0.0_fp
!! !
!! klc = 1
!! do k = kfsmax(nm),kfsmin(nm),-1
!! concin3d(klc) = max(0.0_fp , r0(nm,k,ll))
!! klc=klc+1
!! enddo
!! !
!! ! Solve equilibrium concentration vertical and
!! ! integrate over vertical
!! !
!! call eqtran(siglc ,thicklc ,kmaxlc ,wslc ,ltur , &
!! & frac(nm,l),tsigmol ,dcwwlc ,lundia ,taucr(l) , &
!! & rksr(nm) ,3 ,lsecfl ,spirint ,suspfrac , &
!! & tetacr(l) ,concin3d , &
!! & dzduu(nm) ,dzdvv(nm) ,ubot(nm) ,tauadd ,sus , &
!! & bed ,susw ,bedw ,espir ,wave , &
!! & scour ,ubot_from_com ,camax ,eps , &
!! & iform(l) ,par(1,l) ,max_integers,max_reals,max_strings, &
!! & dll_function(l),dll_handle(l),dll_integers,dll_reals,dll_strings, &
!! & taks ,caks ,taurat(nm,l),sddflc ,rsdqlc , &
!! & kmaxsd ,conc2d ,sbcu(nm,l ),sbcv(nm,l),sbwu(nm,l), &
!! & sbwv(nm,l),sswu(nm,l),sswv(nm,l),tdss ,caks_ss3d , &
!! & aks_ss3d ,ust2(nm) ,tsd ,error )
!! if (error) call d3stop(1, gdp)
!! if (suspfrac) then
!! aks(nm, l) = taks
!! dss(nm, l) = tdss
!! !
!! ! Copy results into arrays
!! !
!! kmxsed(nm, l) = kfsmin(nm)+kmaxlc-kmaxsd
!! !
!! klc=0
!! do k = kfsmax(nm),kfsmin(nm)-1,-1
!! seddif(nm,k,l) = sddflc(klc)
!! klc = klc + 1
!! enddo
!! klc = 1
!! do k = kfsmax(nm),kfsmin(nm),-1
!! rsedeq(nm,k,l) = rsdqlc(klc)
!! klc = klc + 1
!! enddo
!! !
!! ! Source and sink terms for main 3d computation
!! ! note: terms are part explicit, part implicit, see
!! ! thesis of Giles Lesser, May 2000
!! !
!! thick0 = thicklc(kmaxsd) * h0
!! thick1 = thicklc(kmaxsd) * h1
!! call soursin_3d (h1 ,thick0 ,thick1 , &
!! & siglc(kmaxsd) ,thicklc(kmaxsd),r0(nm,kmxsed(nm,l) ,ll) , &
!! & vicmol ,sigmol(ll) ,seddif(nm,kmxsed(nm,l)-1,l), &
!! & rhosol(l) ,caks_ss3d ,ws(nm,kmxsed(nm,l),l) , &
!! & aks_ss3d ,sourse(nm,l) ,sour_im(nm,l) , &
!! & sinkse(nm,l) )
!! ! Impose relatively large vertical diffusion
!! ! coefficients for sediment in layer interfaces from
!! ! bottom of reference cell downwards, to ensure little
!! ! gradient in sed. conc. exists in this area.
!! !
!! difbot = 10.0_fp * ws(nm,kmxsed(nm,l)-1,l) * thick1
!! do k = kfsmin(nm)-1, kmxsed(nm,1)-1
!! seddif(nm, k, l) = difbot
!! enddo
!! endif ! suspfrac
else
!
! kmaxlc = 1
! 2D CASE (Numerical approximation)
!
if (suspfrac) then
!
! Fill local 1dv arrays with fall velocity and
! diffusivity
!
do k2d = 0, kmax2d
!! ws2d(k2d) = ws(nm, kbed, l)
! BEGIN DEBUG
ws2d(k2d) = 1d-4
! END DEBUG
dcww2d(k2d) = 0.0_fp
enddo
!! trsedeq = rsedeq(nm, kbed, l)
else
trsedeq = 0.0_fp
endif
taks = taks0
!
if (lsecfl > 0) then
!! spirint = r0(nm, kbed, lsecfl)
else
spirint = 0.0_fp
endif
!
! Solve equilibrium concentration vertical and
! integrate over vertical; compute bedload
! transport excluding slope effects.
!
! BEGIN DEBUG
if ( time1.gt.1320 .and. abs(nm-0.5*ndxi).lt.1d0 ) then
continue
end if
! END DEBUG
call eqtran(sig2d ,thck2d ,kmax2d ,ws2d ,ltur , &
& frac(nm,l),tsigmol ,dcww2d ,mdia ,taucr(l) , &
& rksr(nm) ,2 ,lsecfl ,spirint ,suspfrac , &
& tetacr(l) ,concin2d , &
& dzdx ,dzdy ,ubot(nm) ,tauadd ,sus , &
& bed ,susw ,bedw ,espir ,wave , &
& scour ,ubot_from_com ,camax ,eps , &
& iform(l) ,par(:,l) ,max_integers,max_reals,max_strings, &
& dll_function(l),dll_handle(l),dll_integers,dll_reals,dll_strings, &
& taks ,caks ,taurat(nm,l),sddf2d ,rsdq2d , &
& kmaxsd ,trsedeq ,sbcx(nm,l) ,sbcy(nm,l),sbwx(nm,l) , &
& sbwy(nm,l),sswx(nm,l),sswy(nm,l),tdss ,caks_ss3d , &
& aks_ss3d ,ust2(nm) ,tsd ,error )
if (error) return
if (suspfrac) then
!! aks (nm, l) = taks
!! dss (nm, l) = tdss
!! rsedeq(nm, kbed, l) = trsedeq
kmxsed(nm, l) = kbed
!
! Galappatti time scale and source and sink terms
!
call soursin_2d(umod(nm) ,ustarc ,h0 ,h1 , &
!DEBUG & ws(nm,kbed,l) ,tsd ,rsedeq(nm,kbed,l),factsd, &
& twsk ,tsd ,trsedeq, 1d0, &
& sourse(nm,l) ,sour_im(nm,l) ,sinkse(nm,l) )
endif ! suspfrac
endif ! kmaxlc = 1
if (suspfrac) then
! caksrho(l, nm) = caks * rhosol(l)
endif
enddo ! next sediment fraction
enddo ! next nm point
!
! Reduce the source and sink terms to avoid large bed level changes
! Note: previous implementation forgot to multiply source/
! sink terms with the thickness for the 2Dh case
!
!! call z_red_soursin(nmmax ,kmax ,thick ,kmxsed , &
!! & lsal ,ltem ,lsed ,lsedtot , &
!! & dps ,s0 ,s1 ,r0 , &
!! & rsedeq ,nst ,dzs1 ,kfsmax , &
!! & kfsmin ,kfs ,gdp )
!
! Fill sutot and svtot
!
do l = 1,lsedtot
!! call dfexchg( sbcu(:,l) ,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( sbwu(:,l) ,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( sswu(:,l) ,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( sbcv(:,l) ,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( sbwv(:,l) ,1, 1, dfloat, nm_pos, gdp)
!! call dfexchg( sswv(:,l) ,1, 1, dfloat, nm_pos, gdp)
if (sedtyp(l)/=SEDTYP_COHESIVE) then
do nm = 1, ndxi
sxtot(nm, l) = sbcx(nm, l) + sbwx(nm, l) + sswx(nm, l)
sytot(nm, l) = sbcy(nm, l) + sbwy(nm, l) + sswy(nm, l)
enddo
endif
enddo
!
! Upwind scheme for bed load and wave driven transport
! Convert sand bed load transport to U and V points using upwind scheme
!
if (bed > 0.0_fp) then
!
! Upwind bed load transport
!
!! call upwbed(sbcu ,sbcv ,sbcuu ,sbcvv ,kfu , &
!! & kfv ,kcs ,kfsed ,lsedtot , &
!! & nmmax ,icx ,icy ,sutot ,svtot , &
!! & gdp )
call fm_upwbed(lsedtot, sbcx, sbcy, sxtot, sytot, e_sbcn, e_sbct)
endif
!
if (bedw>0.0_fp .and. wave) then
!
! Upwind wave-related bed load load transports
!
!! call upwbed(sbwu ,sbwv ,sbwuu ,sbwvv ,kfu , &
!! & kfv ,kcs ,kfsed ,lsedtot , &
!! & nmmax ,icx ,icy ,sutot ,svtot , &
!! & gdp )
endif
!
if (susw>0.0_fp .and. wave) then
!
! Upwind wave-related suspended load transports
!
!! call upwbed(sswu ,sswv ,sswuu ,sswvv ,kfu , &
!! & kfv ,kcs ,kfsed ,lsedtot , &
!! & nmmax ,icx ,icy ,sutot ,svtot , &
!! & gdp )
endif
!
! Bed-slope and sediment availability effects for
! current-related bed load transport
!
if (bed > 0.0_fp) then
!! call adjust_bedload(nmmax ,icx ,icy ,kcs , &
!! & kcu ,kcv ,kfu ,kfv ,lsedtot , &
!! & sbcuu ,sbcvv ,sbuut ,sbvvt ,dzduu , &
!! & dzdvv ,taurat ,frac ,fixfac ,ust2 , &
!! & hu ,hv ,dm ,hidexp ,.true. , &
!! & .true. ,rhowat ,kmax ,dps ,gsqs , &
!! & guu ,gvv ,guv ,gvu ,kbed , &
!! & gdp )
endif
!
! Bed-slope and sediment availability effects for
! wave-related bed load transport
!
if (bedw>0.0_fp .and. wave) then
!! call adjust_bedload(nmmax ,icx ,icy ,kcs , &
!! & kcu ,kcv ,kfu ,kfv ,lsedtot , &
!! & sbwuu ,sbwvv ,sbuut ,sbvvt ,dzduu , &
!! & dzdvv ,taurat ,frac ,fixfac ,ust2 , &
!! & hu ,hv ,dm ,hidexp ,.true. , &
!! & .false. ,rhowat ,kmax ,dps ,gsqs , &
!! & guu ,gvv ,guv ,gvu ,kbed , &
!! & gdp )
endif
!
! Sediment availability effects for
! wave-related suspended load transport
!
if (susw>0.0_fp .and. wave) then
!! call adjust_bedload(nmmax ,icx ,icy ,kcs , &
!! & kcu ,kcv ,kfu ,kfv ,lsedtot , &
!! & sswuu ,sswvv ,sbuut ,sbvvt ,dzduu , &
!! & dzdvv ,taurat ,frac ,fixfac ,ust2 , &
!! & hu ,hv ,dm ,hidexp ,.false. , &
!! & .false. ,rhowat ,kmax ,dps ,gsqs , &
!! & guu ,gvv ,guv ,gvu ,kbed , &
!! & gdp )
endif
!
! Summation of current-related and wave-related transports
!
!do l = 1,lsedtot
! if (sedtyp(l)/=SEDTYP_COHESIVE) then
! do nm = 1, lnx
! sbuu(l,nm) = e_sbcn(nm, l) + e_sbwn(nm, l) + e_sswn(nm, l)
! enddo
! endif
!enddo
!
! Update sourse fluxes due to sand-mud interaction
!
allocate(E(lsedtot), STAT=istat)
do nm = 1, ndxi
if (pmcrit(nm)<0.0_fp) cycle
if (mudfrac(nm)<=0.0_fp .or. mudfrac(nm)>=1.0_fp) cycle
!
! Compute erosion velocities
!
E = 0.0_fp
do l = 1, lsed
!! ll = lstart + l
kmaxsd = kmxsed (nm,l)
!! if (frac(nm,l)>0.0_fp) E(l) = (sourse(nm,l) - sour_im(nm,l)*r0(nm,kmaxsd,ll))/(cdryb(l)*frac(nm,l))
enddo
!
! Recompute erosion velocities
!
call sand_mud(lsed, E, frac(nm,:), mudfrac(nm), sedtyp, pmcrit(nm))
!
! Recompute erosion fluxes
! only explicit part of erosion flux is changed
! (is effectively the same as changing the equilibrium concentration)
!
do l = 1, lsed
!! ll = lstart + l
kmaxsd = kmxsed (nm,l)
!! sourse(nm,l) = frac(nm,l)*cdryb(l)*E(l) + sour_im(nm,l)*r0(nm,kmaxsd,ll)
enddo
enddo
deallocate(E, STAT=istat)
!
! Add implicit part of source term to sinkse
!
! SPvdP: source_im > 0 acts as sink
do l = 1, lsed
do nm = 1, ndxi
sinkse(nm, l) = sinkse(nm, l) + sour_im(nm, l)
enddo
enddo
! BEGIN DEBUG
! call realloc(plotlin,Ndx)
! plotlin = DMISS
! do nm=1,Ndxi
! plotlin(nm) = sinkse(nm,1)
! plotlin(nm) = sour_im(nm,1)
! plotlin(nm) = sourse(nm,1)
! end do
! END DEBUG
!
! Finally fill sour and sink arrays for both sand and silt
! note that sourse/sinkse and sourf/sinkf arrays are required for BOTT3D
!
do l = 1, lsed
!! ll = ISED1 + l - 1
do nm = 1, ndxi
! k = kmxsed(nm,l)
!
! sourse/sinkse/sourf/sinkf are defined per unit volume
! in z-model the sour/sink terms are defined per content of cell
! sedtra%sourse(nm,l)
!! constsour(ll, k) = constsour(ll, k) + sourse(nm,l)
!! constsink(ll, k) = constsink(ll, k) + sinkse(nm,l)
if (iflufflyr>0) then
!! sour(nm, k, ll) = sour(nm, k, ll) + sourf(l, nm) * gsqs(nm) * dzs0(nm,kmxsed(nm,l))
!! sink(nm, k, ll) = sink(nm, k, ll) + sinkf(l, nm) * gsqs(nm) * dzs1(nm,kmxsed(nm,l))
endif
enddo
!! call dfexchg( sour(:,:,l),1, kmax, dfloat, nm_pos, gdp)
!! call dfexchg( sink(:,:,l),1, kmax, dfloat, nm_pos, gdp)
enddo
!! !
!! ! DD-Mapper: copy sbuu and sbvv
!! !
!! nhystp = nxtstp(d3dflow_sediment, gdp)
return
end subroutine fm_erosed
! Interpolate flownode-based vector (sx,sy) to edge-based vector (sxx, syy)
subroutine fm_upwbed(lsedtot, sx, sy, sxtot, sytot, e_sn, e_st)
use m_flowgeom
use m_flow
use unstruc_messages
implicit none
integer, intent(in) :: lsedtot !< number of sediment fractions
double precision, dimension(Ndx,lsedtot), intent(in) :: sx, sy !< cell (flownode)-based quantity
double precision, dimension(Ndx,lsedtot), intent(in) :: sxtot, sytot !< cell (flownode)-based fluxes
double precision, dimension(Lnx,lsedtot), intent(out) :: e_sn, e_st !< edge (flowlink)-based quantity, normal and tangential components
double precision :: sutot1, sutot2
integer :: k1, k2, Lf, l
!if ( laterallyaveragedbedload ) then
! call mess(LEVEL_ERROR, 'upwbed: laterally averaged bedload not supported')
!end if
! internal flowlinks
do Lf=1,Lnxi
! check if flowlink is active and if it connects two active sediment flownodes
if ( hu(Lf).gt.epshu ) then
! find left and right neighboring flownodes
k1 = ln(1,Lf)
k2 = ln(2,Lf)
do l=1,lsedtot
! project the fluxes in flowlink direction
sutot1 = csu(Lf)*sxtot(k1,l) + snu(Lf)*sytot(k1,l)
sutot2 = csu(Lf)*sxtot(k2,l) + snu(Lf)*sytot(k2,l)
! upwind approximation
if ( sutot1.gt.0d0 .and. sutot2.gt.0d0 ) then
e_sn(Lf,l) = csu(Lf)*sx(k1,l) + snu(Lf)*sy(k1,l)
e_st(Lf,l) = -snu(Lf)*sx(k1,l) + csu(Lf)*sy(k1,l)
else if ( sutot1.lt.0d0 .and. sutot2.lt.0d0 ) then
e_sn(Lf,l) = csu(Lf)*sx(k2,l) + snu(Lf)*sy(k2,l)
e_st(Lf,l) = -snu(Lf)*sx(k2,l) + csu(Lf)*sy(k2,l)
else
e_sn(Lf,l) = csu(Lf)*0.5d0*(sx(k1,l)+sx(k2,l)) + snu(Lf)*0.5d0*(sy(k1,l)+sy(k2,l))
e_st(Lf,l) = -snu(Lf)*0.5d0*(sx(k1,l)+sx(k2,l)) + csu(Lf)*0.5d0*(sy(k1,l)+sy(k2,l))
end if
end do
else
do l=1,lsedtot
e_sn(Lf,l) = 0d0
e_st(Lf,l) = 0d0
end do
end if
end do
! boundary flowlinks
do Lf=Lnxi+1,Lnx
if ( hu(Lf).gt.epshu ) then
! find left and right neighboring flownodes
k1 = ln(1,Lf) ! boundary node
k2 = ln(2,Lf) ! internal node
do l=1,lsedtot
e_sn(Lf,l) = csu(Lf)*sx(k2,l) + snu(Lf)*sy(k2,l)
e_st(Lf,l) = -snu(Lf)*sx(k2,l) + csu(Lf)*sy(k2,l)
end do
else
do l=1,lsedtot
e_sn(Lf,l) = 0d0
e_st(Lf,l) = 0d0
end do
end if
end do
return
end subroutine fm_upwbed
subroutine fm_bott3d()
! nmmax ,kmax ,lsed ,lsedtot , &
! & lsal ,ltem ,kfs ,kfu ,kfv , &
! & r1 ,s0 ,kcs , &
! & dps ,gsqs ,guu , &
! & gvv ,s1 ,thick ,dp , &
! & umean ,vmean ,sbuu ,sbvv , &
! & depchg ,ssuu ,ssvv ,nst ,hu , &
! & hv ,aks ,sig ,u1 ,v1 , &
! & sscomp ,kcsbot , &
! & guv ,gvu ,rca ,kcu , &
! & kcv ,icx ,icy ,timhr , &
! & nto ,volum0 ,volum1 ,dt ,gdp )
!----- GPL ---------------------------------------------------------------------
!
! Copyright (C) Stichting Deltares, 2011-2014.
!
! 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: fm_erosed.f90 42642 2015-10-21 11:34:20Z dam_ar $
! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/fm_erosed.f90 $
!!--description-----------------------------------------------------------------
!
! Function: Computes suspended sediment transport correction
! vector for sand sediment fractions
! Computes depth integrated suspended sediment
! transport vector for output to map file
! Computes change in BODSED based on source and sink
! terms calculated in EROSED, and new concentrations.
! Calculates new mixing layer thickness based on
! change in BODSED values
! Calculates new depth values based on changes
! in bottom sediment.
! Includes erosion of dry points and associated
! bathymetry changes
! Method used: Attention: pointer ll for 'standard' FLOW
! arrays is shifted with lstart
!
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!!--declarations----------------------------------------------------------------
use precision
! use sp_buffer
! use flow_tables
! use bedcomposition_module
! use globaldata
! use dfparall
use sediment_basics_module
use m_flowgeom
use m_sediment, only: stmpar, sedtra, mtd
use m_flowtimes, only: dts, tstart_user, time1
! !
implicit none
! !
! type(globdat),target :: gdp
! !
! ! The following list of pointer parameters is used to point inside the gdp structure
! !
! include 'flow_steps_f.inc'
! integer , pointer :: lundia
real(hp) , pointer :: morft
real(fp) , pointer :: morfac
! real(fp) , pointer :: sus
real(fp) , pointer :: bed
real(fp) , pointer :: tmor
! real(fp) , pointer :: thetsd
! real(fp) , pointer :: sedthr
! real(fp) , pointer :: hmaxth
! integer , pointer :: mergehandle
! integer , pointer :: itmor
! type (handletype) , pointer :: bcmfile
! type (bedbndtype) , dimension(:) , pointer :: morbnd
! real(hp) , dimension(:) , pointer :: mergebuf
logical , pointer :: bedupd
! logical , pointer :: cmpupd
! logical , pointer :: neglectentrainment
! logical , pointer :: multi
! logical , pointer :: wind
! logical , pointer :: temp
! logical , pointer :: const
! logical , pointer :: dredge
! logical , pointer :: struct
! logical , pointer :: sedim
! logical , pointer :: scour
! logical , pointer :: snelli
! real(fp), dimension(:) , pointer :: factor
! real(fp) , pointer :: slope
! real(fp), dimension(:) , pointer :: bc_mor_array
real(fp), dimension(:,:) , pointer :: dbodsd
! real(fp), dimension(:) , pointer :: dm
! real(fp), dimension(:) , pointer :: dg
! real(fp), dimension(:,:) , pointer :: fixfac
! real(fp), dimension(:,:) , pointer :: frac
! integer , dimension(:) , pointer :: kfsed
! integer , dimension(:,:) , pointer :: kmxsed
! real(fp), dimension(:) , pointer :: mudfrac
! real(fp), dimension(:,:) , pointer :: sbuuc
! real(fp), dimension(:,:) , pointer :: sbvvc
! real(fp), dimension(:,:) , pointer :: ssuuc
! real(fp), dimension(:,:) , pointer :: ssvvc
! real(fp), dimension(:,:) , pointer :: sucor
! real(fp), dimension(:,:) , pointer :: svcor
real(fp), dimension(:,:) , pointer :: sinkse
real(fp), dimension(:,:) , pointer :: sourse
! integer , pointer :: nmudfrac
! real(fp) , dimension(:) , pointer :: rhosol
real(fp) , dimension(:) , pointer :: cdryb
integer , dimension(:) , pointer :: sedtyp
! integer , pointer :: julday
! integer , pointer :: ntstep
! real(fp), dimension(:,:,:) , pointer :: fluxu
! real(fp), dimension(:,:,:) , pointer :: fluxv
! real(fp), dimension(:) , pointer :: duneheight
! integer , pointer :: iflufflyr
! real(fp), dimension(:,:) , pointer :: mfluff
! real(fp), dimension(:,:) , pointer :: sinkf
! real(fp), dimension(:,:) , pointer :: sourf
!!
!! Local parameters
!!
! integer, parameter :: bedchangemessmax = 50
!!
!! Global variables
!!
! integer , intent(in) :: icx
! integer , intent(in) :: icy
! integer , intent(in) :: kmax ! Description and declaration in esm_alloc_int.f90
! integer , intent(in) :: lsal ! Description and declaration in dimens.igs
! integer , intent(in) :: lsed ! Description and declaration in esm_alloc_int.f90
! integer , intent(in) :: lsedtot! Description and declaration in esm_alloc_int.f90
integer, pointer :: lsedtot
integer, pointer :: lsed
! integer , intent(in) :: ltem ! Description and declaration in dimens.igs
! integer , intent(in) :: nmmax ! Description and declaration in dimens.igs
! integer , intent(in) :: nto ! Number of open boundaries (esm_alloc_int.igs)
! integer , intent(in) :: nst
! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kcs ! Description and declaration in esm_alloc_int.f90
! integer , dimension(gdp%d%nmlb:gdp%d%nmub) :: kcsbot
! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kcu ! Description and declaration in esm_alloc_int.f90
! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kcv ! Description and declaration in esm_alloc_int.f90
! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfs ! Description and declaration in esm_alloc_int.f90
! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfu ! Description and declaration in esm_alloc_int.f90
! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfv ! Description and declaration in esm_alloc_int.f90
! logical , intent(in) :: sscomp
! real(fp) , intent(in) :: dt
! real(fp) :: timhr
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, lsed) , intent(in) :: aks ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: depchg ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(:), pointer :: depchg
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: dp ! Description and declaration in esm_alloc_real.f90
! real(prec), dimension(gdp%d%nmlb:gdp%d%nmub) :: dps ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gsqs ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: guu ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: guv ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gvu ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gvv ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: hu ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: hv ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: s0 ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) :: s1 ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: umean ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: vmean ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: u1 ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: v1 ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: volum0 ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: volum1 ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, *), intent(in) :: r1 ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, lsed) , intent(in) :: rca ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, lsedtot) :: sbuu ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, lsedtot) :: sbvv ! Description and declaration in esm_alloc_real.f90
real(fp), dimension(:,:), pointer :: e_sbcn
real(fp), dimension(:,:), pointer :: e_sbct
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, lsed) :: ssuu ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, lsed) :: ssvv ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(kmax) , intent(in) :: sig ! Description and declaration in esm_alloc_real.f90
! real(fp), dimension(kmax) , intent(in) :: thick ! Description and declaration in esm_alloc_real.f90
!!
!! Local variables
!!
! integer :: i
! integer :: ib
! integer :: icond
! integer :: idir_scalar
! integer :: jb
! integer :: k
! integer :: kvalue
integer :: l
! integer :: li
! integer :: ll
! integer :: lsedbed
! integer :: lstart
! integer :: m
! integer :: n
! integer :: ndm
! integer :: nhystp
integer :: nm
! integer :: nmd
! integer :: nm_pos ! indicating the array to be exchanged has nm index at the 2nd place, e.g., dbodsd(lsedtot,nm)
! integer :: nmu
! integer :: num
! integer :: numu
! integer :: nxmx
! integer :: bedchangemesscount
logical :: bedload
! logical :: from_ndm
! logical :: from_nmd
! logical :: from_nmu
! logical :: from_num
! real(fp) :: aksu
! real(fp) :: apower
! real(fp) :: cavg
! real(fp) :: cavg1
! real(fp) :: cavg2
! real(fp) :: ceavg
! real(fp) :: cumflux
! real(fp) :: dhmax
! real(fp) :: alfa_dist
! real(fp) :: alfa_mag
real(fp) :: dsdnm
! real(fp) :: dv
! real(fp) :: dz
real(fp) :: eroflx
! real(fp) :: fact
! real(fp) :: gsqsmin
! real(fp) :: gsqsinv
! real(fp) :: h1
real(fp) :: dtmor
! real(fp) :: htdif
! real(fp) :: rate
! real(fp) :: r1avg
real(fp) :: sedflx
! real(fp) :: thet
! real(fp) :: thick0
! real(fp) :: thick1
! real(fp) :: totdbodsd
! real(fp) :: totfixfrac
real(fp) :: trndiv
! real(fp) :: z
real(fp) :: flux, sumflux
integer :: ii, LL, Lf
!!
!!! executable statements -------------------------------------------------------
!!
! lundia => gdp%gdinout%lundia
morft => stmpar%morpar%morft
morfac => stmpar%morpar%morfac
! sus => gdp%gdmorpar%sus
! bed => gdp%gdmorpar%bed
bed => stmpar%morpar%bed
tmor => stmpar%morpar%tmor
! thetsd => gdp%gdmorpar%thetsd
! sedthr => gdp%gdmorpar%sedthr
! hmaxth => gdp%gdmorpar%hmaxth
! mergehandle => gdp%gdmorpar%mergehandle
! itmor => gdp%gdmorpar%itmor
! bcmfile => gdp%gdmorpar%bcmfile
! morbnd => gdp%gdmorpar%morbnd
! mergebuf => gdp%gdmorpar%mergebuf
bedupd => stmpar%morpar%bedupd
! cmpupd => gdp%gdmorpar%cmpupd
! neglectentrainment => gdp%gdmorpar%neglectentrainment
! multi => gdp%gdmorpar%multi
! wind => gdp%gdprocs%wind
! temp => gdp%gdprocs%temp
! const => gdp%gdprocs%const
! dredge => gdp%gdprocs%dredge
! struct => gdp%gdprocs%struct
! sedim => gdp%gdprocs%sedim
! snelli => gdp%gdprocs%snelli
! scour => gdp%gdscour%scour
! factor => gdp%gdscour%factor
! slope => gdp%gdscour%slope
! bc_mor_array => gdp%gderosed%bc_mor_array
dbodsd => sedtra%dbodsd
! dm => gdp%gderosed%dm
! dg => gdp%gderosed%dg
! fixfac => gdp%gderosed%fixfac
! frac => gdp%gderosed%frac
! kfsed => gdp%gderosed%kfsed
! kmxsed => gdp%gderosed%kmxsed
! mudfrac => gdp%gderosed%mudfrac
! sbuuc => gdp%gderosed%e_sbnc
! sbvvc => gdp%gderosed%e_sbtc
! ssuuc => gdp%gderosed%e_ssnc
! ssvvc => gdp%gderosed%e_sstc
! sucor => gdp%gderosed%e_scrn
! svcor => gdp%gderosed%e_scrt
sinkse => sedtra%sinkse
sourse => sedtra%sourse
! nmudfrac => gdp%gdsedpar%nmudfrac
! rhosol => gdp%gdsedpar%rhosol
cdryb => stmpar%sedpar%cdryb
! sedtyp => gdp%gdsedpar%sedtyp
sedtyp => stmpar%sedpar%sedtyp
! julday => gdp%gdinttim%julday
! ntstep => gdp%gdinttim%ntstep
! fluxu => gdp%gdflwpar%fluxu
! fluxv => gdp%gdflwpar%fluxv
! duneheight => gdp%gdbedformpar%duneheight
! iflufflyr => gdp%gdmorpar%flufflyr%iflufflyr
! mfluff => gdp%gdmorpar%flufflyr%mfluff
! sinkf => gdp%gdmorpar%flufflyr%sinkf
! sourf => gdp%gdmorpar%flufflyr%sourf
lsed => stmpar%lsedsus
lsedtot => stmpar%lsedtot
e_sbcn => sedtra%e_sbcn
e_sbct => sedtra%e_sbct
depchg => mtd%depchg
! !
! lstart = max(lsal, ltem)
bedload = .false.
dtmor = dts*morfac
! nm_pos = 1
! !
! ! parallel case: exchange arrays in the overlapping regions
! ! other arrays are their overlap data exchanged in erosed.f90
! !
! call dfexchg( fixfac,1, lsedtot, dfloat, nm_pos, gdp)
! call dfexchg( frac, 1, lsedtot, dfloat, nm_pos, gdp)
! if (lsed > 0) then
! call dfexchg( aks, 1, lsed, dfloat, nm_pos, gdp)
! call dfexchg( rca, 1, lsed, dfloat, nm_pos, gdp)
! endif
! !
! ! Calculate suspended sediment transport correction vector (for SAND)
! ! Note: uses GLM velocites, consistent with DIFU
! !
! ! Correct suspended sediment transport rates by estimating the
! ! quantity of suspended sediment transported in the grid cells below
! ! Van Rijn's reference height (aks) and making a vector of this in the
! ! opposite direction to the suspended sediment transport.
! !
! ! ensure suspended sediment correction arrays and suspended sediment
! ! vector arrays are blank
! !
! !
! sucor = 0.0_fp
! svcor = 0.0_fp
! ssuu = 0.0_fp
! ssvv = 0.0_fp
! !
! ! calculate corrections
! !
! if (sus /= 0.0_fp) then
! !
! ! suspension transport correction vector only for 3D
! !
! if (kmax > 1) then
! do l = 1, lsed
! ll = lstart + l
! if (sedtyp(l) == SEDTYP_NONCOHESIVE_SUSPENDED) then
! !
! call dfexchg( fluxu(:,:,ll) ,1, kmax, dfloat, nm_pos, gdp)
! call dfexchg( fluxv(:,:,ll) ,1, kmax, dfloat, nm_pos, gdp)
! do nm = 1, nmmax
! nmu = nm + icx
! num = nm + icy
! !
! ! try new approach - should be smoother
! ! don't worry about direction of the flow
! ! use concentration at velocity point=average of the
! ! two adjacent concentrations
! ! use aks height at velocity point = average of the
! ! two adjacent aks values
! !
! ! note correction vector only computed for velocity
! ! points with active sediment cells on both sides
! !
! ! u direction
! !
! if ((kfu(nm)*kfsed(nm)*kfsed(nmu)) /= 0) then
! cumflux = 0.0_fp
! if (kcs(nmu) == 3 .or. kcs(nmu) == -1) then
! aksu = aks(nm, l)
! elseif (kcs(nm) == 3 .or. kcs(nm) == -1) then
! aksu = aks(nmu, l)
! else
! aksu = (aks(nm, l) + aks(nmu, l))/2.0
! endif
! !
! ! work up through layers integrating transport
! ! below aksu
! !
! do k = kmax, 1, -1
! kvalue = k
! htdif = aksu/hu(nm) - (1.0 + sig(k) - thick(k)/2.0)
! !
! ! if layer containing aksu
! !
! if (htdif <= thick(k)) then
! cumflux = cumflux + fluxu(nm, k, ll)*htdif/thick(k)
! exit
! else
! cumflux = cumflux + fluxu(nm, k, ll)
! endif
! enddo
! cumflux = cumflux / guu(nm)
! !
! ! integration finished
! ! suspended transport correction = opposite of the
! ! transport below aksu
! !
! ! Approximation of the additional transport between bottom of
! ! kmaxsd layer and za has been included in the correction vector
! !
! k = kvalue
! if (k > 1) then
! if (kcs(nmu) == 3 .or. kcs(nmu) == -1) then
! !
! ! correction for domain decomposition:
! !
! ceavg = rca(nm, l)
! elseif (kcs(nm) == 3 .or. kcs(nm) == -1) then
! ceavg = rca(nmu, l)
! else
! ceavg = (rca(nm, l) + rca(nmu, l))/2.0_fp
! endif
! r1avg = (r1(nm, k - 1, ll) + r1(nmu, k - 1, ll))/2.0
! if (ceavg>r1avg*1.1 .and. ceavg>0.05) then
! z = (1.0 + sig(k - 1))*hu(nm)
! apower = log(max(r1avg/ceavg,1.0e-5_fp))/log(z/aksu)
! z = (1.0 + sig(k - 1) - 0.5*thick(k - 1))*hu(nm)
! dz = (thick(k)-htdif)*hu(nm)
! if (apower>-1.05 .and. apower<=-1.0) then
! apower = -1.05
! elseif (apower>=-1.0 .and. apower<-0.95) then
! apower = -0.95
! else
! endif
! apower = max(-10.0_fp , apower)
! cavg1 = (ceavg/(apower+1.0)) * (1.0/aksu)**apower
! cavg2 = z**(apower+1.0) - aksu**(apower+1.0)
! cavg = cavg1 * cavg2 / dz
! cumflux= cumflux - u1(nm,k)*(cavg-r1avg)*dz
! endif
! endif
! sucor(nm,l) = -cumflux
! !
! ! bedload will be reduced in case of sediment transport
! ! over a non-erodible layer (no sediment in bed) in such
! ! a case, the suspended sediment transport vector must
! ! also be reduced.
! !
! if ((sucor(nm,l)>0.0_fp .and. kcs(nm)==1) .or. kcs(nmu)/=1) then
! sucor(nm,l) = sucor(nm,l) * fixfac(nm,l)
! else
! sucor(nm,l) = sucor(nm,l) * fixfac(nmu,l)
! endif
! endif
! !
! ! v direction
! !
! if ((kfv(nm)*kfsed(nm)*kfsed(num)) /= 0) then
! cumflux = 0.0_fp
! if (kcs(num) == 3 .or. kcs(num) == -1) then
! aksu = aks(nm, l)
! elseif (kcs(nm) == 3 .or. kcs(nm) == -1) then
! aksu = aks(num, l)
! else
! aksu = (aks(nm, l)+aks(num, l)) / 2.0_fp
! endif
! !
! ! work up through layers integrating transport
! ! below aksu
! !
! do k = kmax, 1, -1
! kvalue = k
! htdif = aksu/hv(nm) - (1.0 + sig(k) - thick(k)/2.0)
! !
! ! if layer containing aksu
! !
! if (htdif <= thick(k)) then
! cumflux = cumflux + fluxv(nm,k,ll)*htdif/thick(k)
! exit
! else
! cumflux = cumflux + fluxv(nm,k,ll)
! endif
! enddo
! cumflux = cumflux / gvv(nm)
! !
! ! integration finished
! ! suspended transport correction = opposite of the
! ! transport below aksu
! !
! ! Approximation of the additional transport between bottom of
! ! kmaxsd layer and za has been included in the correction vector
! !
! k = kvalue
! if (k > 1) then
! if (kcs(num) == 3 .or. kcs(num) == -1) then
! !
! ! correction for domain decomposition:
! !
! ceavg = rca(nm,l)
! elseif (kcs(nm) == 3 .or. kcs(nm) == -1) then
! ceavg = rca(num,l)
! else
! ceavg = (rca(nm,l)+rca(num,l)) / 2.0
! endif
! r1avg = (r1(nm,k-1,ll)+r1(num,k-1,ll)) / 2.0
! if (ceavg>r1avg*1.1 .and. ceavg>0.05) then
! z = (1.0+sig(k-1)) * hv(nm)
! apower = log(max(r1avg/ceavg,1.0e-5_fp))/log(z/aksu)
! z = (1.0+sig(k-1)-0.5*thick(k-1)) * hv(nm)
! dz = (thick(k)-htdif) * hv(nm)
! if (apower>-1.05 .and. apower<=-1.0) then
! apower = -1.05
! elseif (apower>=-1.0 .and. apower<-0.95) then
! apower = -0.95
! else
! endif
! apower = max(-10.0_fp , apower)
! cavg1 = (ceavg/(apower+1.0)) * (1/aksu)**apower
! cavg2 = z**(apower+1.0) - aksu**(apower+1.0)
! cavg = cavg1 * cavg2 / dz
! cumflux= cumflux - v1(nm,k)*(cavg-r1avg)*dz
! endif
! endif
! svcor(nm, l) = -cumflux
! !
! ! bedload will be reduced in case of sediment transport
! ! over a non-erodible layer (no sediment in bed) in such
! ! a case, the suspended sediment transport vector must
! ! also be reduced.
! !
! if ((svcor(nm,l) > 0.0_fp .and. kcs(nm)==1) .or. kcs(num)/=1) then
! svcor(nm, l) = svcor(nm, l)*fixfac(nm, l)
! else
! svcor(nm, l) = svcor(nm, l)*fixfac(num, l)
! endif
! endif
! enddo ! nm
! endif ! sedtyp = SEDTYP_NONCOHESIVE_SUSPENDED
! enddo ! l
! endif ! kmax>1
! !
! if (lsed > 0) then
! call dfexchg( sucor, 1, lsed, dfloat, nm_pos, gdp)
! call dfexchg( svcor, 1, lsed, dfloat, nm_pos, gdp)
! endif
! !
! ! Calculate suspended sediment transport vector components for
! ! output
! ! Note: uses DIFU fluxes
! ! if suspended sediment vector is required this half timestep
! ! note, will be required if nst.ge.itmor for cumulative
! ! transports
! !
! if (sscomp .or. nst>=itmor) then
! do l = 1, lsed
! ll = lstart + l
! do nm = 1, nmmax
! nmu = nm + icx
! num = nm + icy
! !
! ! u component
! !
! if (kfu(nm) == 1 .and. kcu(nm)/=-1) then
! cumflux = 0.0_fp
! do k = 1, kmax
! cumflux = cumflux + fluxu(nm, k, ll)
! enddo
! !
! ! total suspended transport
! !
! ssuu(nm, l) = cumflux/guu(nm) + sucor(nm, l)
! endif
! !
! ! v component
! !
! if (kfv(nm) == 1 .and. kcv(nm)/=-1) then
! cumflux = 0.0_fp
! do k = 1, kmax
! cumflux = cumflux + fluxv(nm, k, ll)
! enddo
! !
! ! total suspended transport
! !
! ssvv(nm, l) = cumflux/gvv(nm) + svcor(nm, l)
! endif
! enddo ! nm
! enddo ! l
! endif ! sscomp .or. nst>=itmor
! endif ! sus /= 0.0
! !
! if (lsed > 0) then
! call dfexchg( ssuu, 1, lsed, dfloat, nm_pos, gdp)
! call dfexchg( ssvv, 1, lsed, dfloat, nm_pos, gdp)
! endif
!
! if morphological computations have started
!
! if (nst >= itmor)
if (time1 > tstart_user + tmor * 60._fp) then ! tmor in minutes since start of computations, time1 in seconds since reference date
! !
! ! Increment morphological time
! ! Note: dtmor in seconds, morft in days!
! !
morft = morft + dtmor/86400.0_fp
! !
! ! Bed boundary conditions: transport condition
! !
! do jb = 1, nto
! icond = morbnd(jb)%icond
! if (icond == 4 .or. icond == 5) then
! !
! ! Open boundary with transport boundary condition:
! ! Get data from table file
! !
! call flw_gettabledata(bcmfile , morbnd(jb)%ibcmt(1) , &
! & morbnd(jb)%ibcmt(2) , morbnd(jb)%ibcmt(3) , &
! & morbnd(jb)%ibcmt(4) , bc_mor_array , &
! & timhr ,julday , gdp )
! !
! ! Prepare loop over boundary points
! !
! do ib = 1, morbnd(jb)%npnt
! alfa_dist = morbnd(jb)%alfa_dist(ib)
! alfa_mag = morbnd(jb)%alfa_mag(ib)
! idir_scalar = morbnd(jb)%idir(ib)
! nm = morbnd(jb)%nm(ib)
! nxmx = morbnd(jb)%nxmx(ib)
! !
! nmu = nm + icx
! num = nm + icy
! nmd = nm - icx
! ndm = nm - icy
! !
! ! If the computed transport is directed outward, do not
! ! impose the transport rate (at outflow boundaries the
! ! "free bed level boundary" condition is imposed. This
! ! check is carried out for each individual boundary point.
! !
! ! Detect the case based on the value of nxmx.
! !
! if (nxmx == nmu) then
! if (umean(nm)<0.0_fp) cycle
! elseif (nxmx == nmd) then
! if (umean(nmd)>0.0_fp) cycle
! elseif (nxmx == num) then
! if (vmean(nm)<0.0_fp) cycle
! elseif (nxmx == ndm) then
! if (vmean(ndm)>0.0_fp) cycle
! endif
! !
! ! The velocity/transport points to the left and top are part
! ! of this cell. nxmx contains by default the index of the
! ! neighbouring grid cell, so that has to be corrected. This
! ! correction is only carried out locally since we need the
! ! unchanged nxmx value further down for the bed level updating
! !
! if (nxmx == nmu .or. nxmx == num) nxmx = nm
! !
! li = 0
! lsedbed = lsedtot - nmudfrac
! do l = 1, lsedtot
! !
! ! bed load transport always zero for mud fractions
! !
! if (sedtyp(l) == SEDTYP_COHESIVE) cycle
! li = li + 1
! !
! if (morbnd(jb)%ibcmt(3) == lsedbed) then
! rate = bc_mor_array(li)
! elseif (morbnd(jb)%ibcmt(3) == 2*lsedbed) then
! rate = bc_mor_array(li) + &
! & alfa_dist * (bc_mor_array(li+lsedbed)-bc_mor_array(li))
! endif
! rate = alfa_mag * rate
! !
! if (icond == 4) then
! !
! ! transport including pores
! !
! rate = rate*cdryb(l)
! else
! !
! ! transport excluding pores
! !
! rate = rate*rhosol(l)
! endif
! !
! ! impose boundary condition
! !
! if (idir_scalar == 1) then
! sbuu(nxmx, l) = rate
! else
! sbvv(nxmx, l) = rate
! endif
! enddo ! l (sediment fraction)
! enddo ! ib (boundary point)
! endif ! icond = 4 or 5 (boundary with transport condition)
! enddo ! jb (open boundary)
! !
! call dfexchg(sbuu, 1, lsedtot, dfloat, nm_pos, gdp)
! call dfexchg(sbvv, 1, lsedtot, dfloat, nm_pos, gdp)
!
! Update quantity of bottom sediment
!
dbodsd = 0.0_fp
!
! compute change in bodsed (dbodsd)
!
! bedchangemesscount = 0
do l = 1, lsedtot
bedload = sedtyp(l)==SEDTYP_NONCOHESIVE_TOTALLOAD
! ll = lstart + l
! do nm = 1, nmmax
do nm=1,Ndxi
! !
! ! note: do not update bottom sediment at open boundary pnts
! !
! if (kcs(nm)*kfs(nm) /= 1) cycle
! !
! nmu = nm + icx
! num = nm + icy
! nmd = nm - icx
! ndm = nm - icy
trndiv = 0.0_fp
sedflx = 0.0_fp
eroflx = 0.0_fp
! gsqsinv = 1.0_fp/gsqs(nm)
! if (sus/=0.0_fp .and. .not. bedload) then
! if (neglectentrainment) then
! !
! ! mass balance based on fluxes: entrainment and deposition
! ! does not lead to erosion/sedimentation.
! !
! if (snelli) then
! !
! ! Only cross-shore component
! !
! trndiv = trndiv + gsqsinv &
! & *( ssuu(nmd, l)*guu(nmd) - ssuu(nm, l)*guu(nm))
! else
! trndiv = trndiv + gsqsinv &
! & *( ssuu(nmd, l)*guu(nmd) - ssuu(nm, l)*guu(nm) &
! & + ssvv(ndm, l)*gvv(ndm) - ssvv(nm, l)*gvv(nm))
! endif
! else
! !
! ! mass balance includes entrainment and deposition
! !
! if (sedtyp(l) == SEDTYP_NONCOHESIVE_SUSPENDED) then
! !
! ! l runs from 1 to lsedtot, kmxsed is defined for 1:lsed
! ! The first lsed fractions are the suspended fractions,
! ! so this goes right
! !
! k = kmxsed(nm, l)
! else
! k = kmax
! endif
! thick0 = volum0(nm,k) * gsqsinv
! thick1 = volum1(nm,k) * gsqsinv
! sedflx = sinkse(nm,l) * r1(nm,k,ll) * thick1
! eroflx = sourse(nm,l) * thick0
! !
! ! Update fluff layer
! !
! if (iflufflyr>0) then
! mfluff(l, nm) = mfluff(l, nm) + &
! & dt*( sinkf(l, nm)*r1(nm, k, ll)*thick1 &
! & - sourf(l, nm) *thick0 )
! endif
! !
! ! add suspended transport correction vector
! !
! if (snelli) then
! !
! ! Only cross-shore component
! !
! trndiv = trndiv + gsqsinv &
! & *( sucor(nmd,l)*guu(nmd)-sucor(nm,l)*guu(nm))
! else
! trndiv = trndiv + gsqsinv &
! & *( sucor(nmd,l)*guu(nmd)-sucor(nm,l)*guu(nm) &
! & +svcor(ndm,l)*gvv(ndm)-svcor(nm,l)*gvv(nm))
! endif
! endif
! endif
if (bed /= 0.0_fp) then
! if (snelli) then
! !
! ! Only cross-shore component
! !
! trndiv = trndiv + gsqsinv &
! & *( sbuu(nmd,l)*guu(nmd)-sbuu(nm,l)*guu(nm))
! else
! trndiv = trndiv + gsqsinv &
! & *( sbuu(nmd,l)*guu(nmd)-sbuu(nm,l)*guu(nm) &
! & +sbvv(ndm,l)*gvv(ndm)-sbvv(nm,l)*gvv(nm))
sumflux = 0d0
do ii=1,nd(nm)%lnx
LL = nd(nm)%ln(ii)
Lf = iabs(LL)
flux = e_sbcn(Lf,l)*wu(Lf)
if ( LL.gt.0 ) then ! inward
sumflux = sumflux + flux
else ! outward
sumflux = sumflux - flux
end if
end do
trndiv = trndiv + sumflux * bai(nm)
! endif
endif
! !
dsdnm = (trndiv+sedflx-eroflx) * dtmor
! !
! ! Warn if bottom changes are very large,
! ! depth change NOT LIMITED
! !
! dhmax = 0.05_fp
! h1 = max(0.01_fp, s1(nm) + real(dps(nm),fp))
! if (abs(dsdnm) > dhmax*h1*cdryb(1) .and. bedupd) then
! !
! ! Only write bed change warning when bed updating is true
! ! (otherwise no problem)
! ! Limit the number of messages with bedchangemessmax
! ! (otherwise tri-diag will grow very fast)
! !
! bedchangemesscount = bedchangemesscount + 1
! if (bedchangemesscount <= bedchangemessmax) then
! call nm_to_n_and_m(nm, n, m, gdp)
! write (lundia, '(a,f5.1,a,i0,a,i0,a,i0,a)') &
! & '*** WARNING Bed change exceeds ' , dhmax*100.0_fp, ' % of waterdepth after ', ntstep, &
! & ' timesteps, location (m,n) = (', m,',',n,')'
! endif
! endif
! !
! ! Update dbodsd value at nm
! !
dbodsd(l, nm) = dbodsd(l, nm) + dsdnm
! !
! call updwaqflxsed(nst, nm, l, trndiv, sedflx, eroflx, gdp)
enddo ! nm
enddo ! l
! if (bedchangemesscount > bedchangemessmax) then
! write (lundia,'(12x,a,i0,a)') 'Bed change messages skipped (more than ',bedchangemessmax,')'
! write (lundia,'(12x,2(a,i0))') 'Total number of Bed change messages for timestep ',ntstep,' : ',bedchangemesscount
! endif
! !
! call fluff_burial(gdp%gdmorpar%flufflyr, dbodsd, lsed, lsedtot, gdp%d%nmlb, gdp%d%nmub, dt, morfac)
! !
! nm_pos = 2
! call dfexchg(dbodsd, 1, lsedtot, dfloat, nm_pos, gdp)
! nm_pos = 1
! !
! ! Re-distribute erosion near dry and shallow points to allow erosion
! ! of dry banks
! !
! do nm = 1, nmmax
! !
! ! If this is a cell in which sediment processes are active then ...
! !
! if (kcs(nm)*kfs(nm)*kfsed(nm) /= 1) cycle
! !
! nmu = nm + icx
! num = nm + icy
! nmd = nm - icx
! ndm = nm - icy
! totdbodsd = 0.0_fp
! do l = 1, lsedtot
! totdbodsd = totdbodsd + dbodsd(l, nm)
! enddo
! !
! ! If this is a cell in erosion is occuring (accretion is not
! ! distributed to dry points) then...
! !
! if (totdbodsd < 0.0_fp) then
! !
! ! Note: contrary to the previous implementation, this new
! ! implementation erodes the sediment from nm and
! ! re-distributes the eroded volume based on the composition
! ! of the neighbouring cells, replenishing the sediment volume
! ! at grid point nm with sediment of a different composition
! ! than that what was eroded. This new implementation is mass
! ! conserving per fraction. Furthermore, re-distribution takes
! ! place only in case of net TOTAL erosion, i.e. not of
! ! individual fractions.
! !
! gsqsmin = gsqs(nm)
! totfixfrac = 0.0_fp
! !
! from_ndm = kfsed(ndm)==0 .and. kcs(ndm) /= 0 .and. kcs(ndm)<3 .and. kcv(ndm)==1 .and. dps(ndm) 1.0e-7_fp) then
! !
! ! Compute local re-distribution factor THET
! !
! if (hmaxth > sedthr) then
! h1 = real(dps(nm),fp) + s1(nm)
! thet = (h1 - sedthr)/(hmaxth - sedthr)*thetsd
! thet = min(thet, thetsd)
! else
! thet = thetsd
! endif
! !
! ! Combine some constant factors in variable THET
! ! Note: TOTDBODSD<0.0 and thus THET>0.0 !
! !
! thet = -gsqsmin * totdbodsd * thet / totfixfrac
! !
! do l = 1, lsedtot
! !
! ! update dbodsd values in this cell and surrounding cells
! ! adjust bedload transport rates to include this erosion
! ! process.
! !
! if (from_ndm) then
! dv = thet*fixfac(ndm, l)*frac(ndm, l)
! dbodsd(l, ndm) = dbodsd(l, ndm) - dv/gsqs(ndm)
! dbodsd(l, nm ) = dbodsd(l, nm ) + dv/gsqs(nm)
! sbvv(ndm, l) = sbvv(ndm, l) + dv/(dtmor*gvv(ndm))
! endif
! !
! if (from_nmd) then
! dv = thet*fixfac(nmd, l)*frac(nmd, l)
! dbodsd(l, nmd) = dbodsd(l, nmd) - dv/gsqs(nmd)
! dbodsd(l, nm ) = dbodsd(l, nm ) + dv/gsqs(nm)
! sbuu(nmd, l) = sbuu(nmd, l) + dv/(dtmor*guu(nmd))
! endif
! !
! if (from_nmu) then
! dv = thet*fixfac(nmu, l)*frac(nmu, l)
! dbodsd(l, nmu) = dbodsd(l, nmu) - dv/gsqs(nmu)
! dbodsd(l, nm ) = dbodsd(l, nm ) + dv/gsqs(nm)
! sbuu(nm, l) = sbuu(nm, l) - dv/(dtmor*guu(nm))
! endif
! !
! if (from_num) then
! dv = thet*fixfac(num, l)*frac(num, l)
! dbodsd(l, num) = dbodsd(l, num) - dv/gsqs(num)
! dbodsd(l, nm ) = dbodsd(l, nm ) + dv/gsqs(nm)
! sbvv(nm, l) = sbvv(nm, l) - dv/(dtmor*gvv(nm))
! endif
! enddo ! l
! endif ! totfixfrac > 1.0e-7
! endif ! totdbodsd < 0.0
! enddo ! nm
! !
! nm_pos = 2
! call dfexchg(dbodsd, 1, lsedtot, dfloat, nm_pos, gdp)
! nm_pos = 1
! !
! ! Modifications for running parallel conditions
! !
! !
! if (multi) then
! i = 0
! do l = 1, lsedtot
! do nm = 1, nmmax
! i = i + 1
! mergebuf(i) = real(dbodsd(l, nm),hp)
! enddo
! enddo
! call putarray (mergehandle,mergebuf(1:nmmax*lsedtot),nmmax*lsedtot)
! call getarray (mergehandle,mergebuf(1:nmmax*lsedtot),nmmax*lsedtot)
! i = 0
! do l = 1, lsedtot
! do nm = 1, nmmax
! i = i + 1
! dbodsd(l, nm) = real(mergebuf(i),fp)
! enddo
! enddo
! endif
! !
! ! Add transports to cumulative transports
! !
! do l = 1, lsedtot
! do nm = 1, nmmax
! sbuuc(nm, l) = sbuuc(nm, l) + sbuu(nm, l) * dtmor
! sbvvc(nm, l) = sbvvc(nm, l) + sbvv(nm, l) * dtmor
! enddo
! enddo
! do l = 1, lsed
! do nm = 1, nmmax
! ssuuc(nm, l) = ssuuc(nm, l) + ssuu(nm, l) * dtmor
! ssvvc(nm, l) = ssvvc(nm, l) + ssvv(nm, l) * dtmor
! enddo
! enddo
! !
! call dfexchg(sbuuc, 1, lsedtot, dfloat, nm_pos, gdp)
! call dfexchg(sbvvc, 1, lsedtot, dfloat, nm_pos, gdp)
! if (lsed > 0) then
! call dfexchg(ssuuc, 1, lsed, dfloat, nm_pos, gdp)
! call dfexchg(ssvvc, 1, lsed, dfloat, nm_pos, gdp)
! endif
! !
! ! Apply erosion and sedimentation to bookkeeping system
! !
! if (cmpupd) then
! !
! ! Determine new thickness of transport layer
! !
! call compthick(dps ,s1 , &
! & nmmax ,gdp )
! !
! ! Update layers and obtain the depth change
! !
! if (updmorlyr(gdp%gdmorlyr, dbodsd, depchg, gdp%messages) /= 0) then
! call writemessages(gdp%messages, lundia)
! call d3stop(1, gdp)
! else
! call writemessages(gdp%messages, lundia)
! endif
! call lyrdiffusion(gdp%gdmorlyr, dtmor)
! !
! ! Apply composition boundary conditions
! !
! call bndmorlyr(lsedtot ,timhr , &
! & nto ,bc_mor_array , &
! & gdp )
! else
! !
! ! Compute bed level changes without actually updating the bed composition
! !
depchg = 0.0_fp
! do nm = 1, nmmax
do nm = 1, Ndxi
! if (kcs(nm)/=0 .and. kcs(nm)<=2) then
do l = 1, lsedtot
depchg(nm) = depchg(nm) + dbodsd(l, nm)/cdryb(l)
enddo
! endif
enddo
! endif
! !
! call dfexchg(depchg, 1, 1, dfloat, nm_pos, gdp)
! !
! ! Bed boundary conditions
! !
! do jb = 1, nto
! icond = morbnd(jb)%icond
! !
! ! In case of an open boundary with bed level condition
! ! described by time series: get data from table file
! !
! if (icond == 2 .or. icond == 3) then
! call flw_gettabledata(bcmfile , morbnd(jb)%ibcmt(1) , &
! & morbnd(jb)%ibcmt(2) , morbnd(jb)%ibcmt(3) , &
! & morbnd(jb)%ibcmt(4) , bc_mor_array , &
! & timhr ,julday , gdp )
! endif
! !
! ! Prepare loop over boundary points
! !
! do ib = 1, morbnd(jb)%npnt
! alfa_dist = morbnd(jb)%alfa_dist(ib)
! alfa_mag = morbnd(jb)%alfa_mag(ib)**2
! idir_scalar = morbnd(jb)%idir(ib)
! nm = morbnd(jb)%nm(ib)
! nxmx = morbnd(jb)%nxmx(ib)
! !
! nmu = nm + icx
! num = nm + icy
! nmd = nm - icx
! ndm = nm - icy
! !
! ! Bed change in open boundary point
! ! Any boundary condition is changed into a "free bed level
! ! boundary" if the computed transport is directed outward.
! !
! ! Detect the case based on the value of nxmx. In case of a
! ! diagonal water level boundary, there will be two separate
! ! entries in the morbnd structure. The sum of alfa_mag(ib)**2
! ! will be equal to 1.
! !
! if (nxmx == nmu) then
! if (umean(nm)<0.0) icond = 0
! elseif (nxmx == nmd) then
! if (umean(nmd)>0.0) icond = 0
! elseif (nxmx == num) then
! if (vmean(nm)<0.0) icond = 0
! elseif (nxmx == ndm) then
! if (vmean(ndm)>0.0) icond = 0
! endif
! !
! select case(icond)
! case (0,4,5)
! !
! ! outflow or free boundary (0)
! ! or prescribed transport with pores (4)
! ! or prescribed transport without pores (5)
! !
! depchg(nm) = depchg(nm) + depchg(nxmx) * alfa_mag
! case (1)
! !
! ! fixed bed level: no update
! !
! ! depchg(nm) = depchg(nm) + 0.0 * alfa_mag
! case (2)
! !
! ! prescribed depth
! ! temporarily store "bed levels" in variable "rate"
! !
! if (morbnd(jb)%ibcmt(3) == 1) then
! rate = bc_mor_array(1)
! elseif (morbnd(jb)%ibcmt(3) == 2) then
! rate = bc_mor_array(1) + &
! & alfa_dist * (bc_mor_array(2)-bc_mor_array(1))
! endif
! !
! depchg(nm) = depchg(nm) + (real(dps(nm),fp)-rate) * alfa_mag
! case (3)
! !
! ! prescribed depth change rate
! !
! if (morbnd(jb)%ibcmt(3) == 1) then
! rate = bc_mor_array(1)
! elseif (morbnd(jb)%ibcmt(3) == 2) then
! rate = bc_mor_array(1) + &
! & alfa_dist * (bc_mor_array(2)-bc_mor_array(1))
! endif
! !
! depchg(nm) = depchg(nm) - rate * alfa_mag * dtmor
! end select
! enddo ! ib (boundary point)
! enddo ! jb (open boundary)
! else
! !
! ! if morphological computations haven't started yet
! !
! do nm = 1, nmmax
! depchg(nm) = 0.0_fp
! enddo
endif ! nst >= itmor
! !
! call dfexchg(depchg, 1, 1, dfloat, nm_pos, gdp)
! !
! ! Update bottom elevations
! !
if (bedupd) then
! !
! ! note: dps and dp are positive downwards.
! !
! do nm = 1, nmmax
do nm = 1, Ndxi
! !
! ! note: if kcs(nm)=0 then depchg(nm)=0.0
! ! should change to following test because depchg may be small
! ! due to truncation errors
! !
if (abs(depchg(nm)) > 0.0_fp) then
! dps(nm) = dps(nm) - real(depchg(nm),prec)
bl(nm) = bl(nm) + depchg(nm)
endif
enddo
! if (scour) then
! !
! ! -Check bottom slopes and apply an avalance effect if needed
! ! -Depths at waterlevel points (dps) will be updated,
! ! to be used for dpu and dpv
! ! -Depth changes will be added to depchg,to be used for dp
! !
! call avalan(dps ,depchg ,gvu ,guv , &
! & icx ,icy ,gsqs ,kcs ,gdp )
! endif
! do nm = 1, nmmax
! !
! ! note: if kcs(nm)=0 then depchg(nm)=0.0
! ! should change to following test because depchg may be small
! ! due to truncation errors
! !
! if (abs(depchg(nm)) >= 0.0) then
! s1(nm) = max(s1(nm), -real(dps(nm),fp))
! s0(nm) = max(s0(nm), -real(dps(nm),fp))
! !
! ! if dry cells are eroded then bring water level down to
! ! bed or maximum water level in surrounding wet cells
! ! (whichever is higher)
! !
! if (kfs(nm) == 0) then
! s1(nm) = s1(nm) + depchg(nm)
! s0(nm) = s0(nm) + depchg(nm)
! endif
! endif
! !
! ! set flag for updating dp points below (note does not = 2 at
! ! open boundaries)
! !
! if (kcs(nm) == 0) then
! kcsbot(nm) = 0
! else
! kcsbot(nm) = 1
! endif
! enddo
! !
! call dfexchg( kcsbot, 1, 1, dfint, nm_pos, gdp)
! !
! ! Dredging and Dumping
! !
! if (dredge) then
! call dredgedump(dbodsd ,cdryb ,nst ,timhr ,morft , &
! & gdp )
! endif
endif
! ! -----------------------------------------------------------
! ! DD_mapper: copy dps and depchg at zeta points
! ! -----------------------------------------------------------
! nhystp = nxtstp(d3dflow_bottom3d, gdp)
! if (bedupd) then
! !
! ! CALDPU is called after BOTT3D in TRISOL when BEDUPD = TRUE
! ! instead of updating dpu/dpv here
! !
! ! Update dp points
! !
! do nm = 1, nmmax
! nmu = nm + icx
! num = nm + icy
! numu = num + icx
! fact = kcsbot(nm) *gsqs(nm) + kcsbot(num) *gsqs(num) &
! & + kcsbot(nmu)*gsqs(nmu) + kcsbot(numu)*gsqs(numu)
! if (fact > 0.0_fp) then
! dp(nm) = dp(nm) - ( depchg(nm) *gsqs(nm) + depchg(num) *gsqs(num) &
! & + depchg(nmu)*gsqs(nmu) + depchg(numu)*gsqs(numu))/fact
! endif
! enddo
! call dfexchg(dp, 1, 1, dfloat, nm_pos, gdp)
! endif
end subroutine fm_bott3d