Skip to content

Merge sedmor branch: updates to the transport formulas #814

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
hrajagers wants to merge 26 commits into
base: main
Choose a base branch
from
Open

Merge sedmor branch: updates to the transport formulas #814

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
26 commits
Select commit Hold shift + click to select a range
902ccaf
UNST-9434
jreyns Nov 19, 2025
414cdd2
Merge remote-tracking branch 'origin/main' into fm/feature/UNST-9434_…
jreyns Nov 19, 2025
9c983f0
UNST-9434
jreyns Nov 19, 2025
054d6fc
UNST-9434
jreyns Nov 19, 2025
773a0c4
UNST-9434
jreyns Nov 20, 2025
613d30a
Merge remote-tracking branch 'origin/main' into fm/feature/UNST-9434_…
jreyns Nov 20, 2025
106a563
Update src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_ke…
jreyns Dec 22, 2025
0c50bcd
Update src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_da…
jreyns Dec 22, 2025
c108921
Merge remote-tracking branch 'remotes/origin/main' into fm/feature/UN…
hrajagers Jan 28, 2026
c79a179
Fortran styler
hrajagers Jan 28, 2026
7c26555
minor changes adding ".0" and some comment updates
hrajagers Jan 29, 2026
faffc99
review comments part 1
jreyns Feb 4, 2026
eabaed8
UNST-9434 review comments part 2
jreyns Feb 11, 2026
9641a7b
remove superfluous vars
jreyns Feb 19, 2026
e976feb
Merge branch 'fm/feature/UNST-9434_merge_branch_sedmor' of https://gi…
jreyns Feb 19, 2026
43098dc
Merge branch 'main' of https://github.com/Deltares/Delft3d into fm/fe…
jreyns Feb 19, 2026
e0b6bd3
push review comment updates
hrajagers Mar 25, 2026
e0177a1
Merge remote-tracking branch 'origin/main' into fm/feature/UNST-9434_…
jreyns Mar 26, 2026
36be299
UNST-9434
jreyns Mar 26, 2026
66c1ccf
UNST-9434
jreyns Mar 26, 2026
0d74a47
Merge branch 'main' into fm/feature/UNST-9434_merge_branch_sedmor
hrajagers Apr 15, 2026
7e1d761
Fix merge error: subroutine determine_linkbased_cumblchg() was lost
hrajagers Apr 15, 2026
702fda8
temporarily undo all changes outside src/utils_gpl/morphology
hrajagers Apr 16, 2026
d43da29
consistently use fp arguments to build both for fp=dp and fp=sp.
hrajagers Apr 16, 2026
56af017
updates based on review
hrajagers Apr 17, 2026
dd0fe72
Merge branch 'main' into all/feature/UNST-9434_sedmor__traform
hrajagers Apr 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 12 additions & 8 deletions src/utils_gpl/morphology/packages/morphology_io/src/rdtrafrm.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1464,13 +1464,13 @@ subroutine traparams(iform ,name ,nparreq ,nparopt ,parkeyw , &
pardef(2) = 0.3_fp
elseif (iform == 19) then
name = 'Van Thiel / Van Rijn (2008)'
nparopt = 13
nparopt = 14
parkeyw(1) = 'facua'
pardef(1) = 0.1_fp
pardef(1) = 0.0_fp
parkeyw(2) = 'facAs'
pardef(2) = 0.1_fp
pardef(2) = 0.2_fp
parkeyw(3) = 'facSk'
pardef(3) = 0.1_fp
pardef(3) = 0.15_fp
parkeyw(4) = 'waveform'
pardef(4) = 2.0_fp ! 1=ruessink, 2=van thiel
parkeyw(5) = 'sws'
Expand All @@ -1491,16 +1491,18 @@ subroutine traparams(iform ,name ,nparreq ,nparopt ,parkeyw , &
pardef(12) = 30.0_fp
parkeyw(13) = 'cmax'
pardef(13) = 0.1_fp
parkeyw(14) = 'alfad50'
pardef(14) = 0.0_fp

elseif (iform == 20) then
name = 'Soulsby / Van Rijn, XBeach flavour'
nparopt = 14
nparopt = 15
parkeyw(1) = 'facua'
pardef(1) = 0.1_fp
pardef(1) = 0.0_fp
parkeyw(2) = 'facAs'
pardef(2) = 0.1_fp
pardef(2) = 0.2_fp
parkeyw(3) = 'facSk'
pardef(3) = 0.1_fp
pardef(3) = 0.15_fp
parkeyw(4) = 'waveform'
pardef(4) = 2.0_fp ! 1=ruessink, 2=van thiel
parkeyw(5) = 'sws'
Expand All @@ -1523,6 +1525,8 @@ subroutine traparams(iform ,name ,nparreq ,nparopt ,parkeyw , &
pardef(13) = 0.1_fp
parkeyw(14) = 'z0'
pardef(14) = 0.006_fp
parkeyw(15) = 'alfad50'
pardef(15) = 0.0_fp
elseif (iform == 21) then
if (name == ' ') name = 'External subroutine'
nparreq = 0
Expand Down
2 changes: 1 addition & 1 deletion src/utils_gpl/morphology/packages/morphology_kernel/src/bedbc2004.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -293,7 +293,7 @@ subroutine bedbc2004(tp ,rhowat , &
uon = max(1.0e-5_fp , uon)
uoff = max(1.0e-5_fp , uoff)
!
uwbih = (0.5_fp*uon**3.0_fp + 0.5_fp*uoff**3.0_fp)**(1.0_fp/3.0_fp) ! Representative peak orbital velocity
uwbih = (0.5_fp*uon**3.0_fp + 0.5_fp*uoff**3.0_fp)**(1.0_fp/3.0_fp) ! Representative peak orbital velocity magnitude

else if (wform==2) then
! Modification by Marcio Boechat Albernaz
Expand Down
10 changes: 7 additions & 3 deletions src/utils_gpl/morphology/packages/morphology_kernel/src/bedtr2004.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -394,14 +394,18 @@ subroutine bedtr2004(u2dh ,d50 ,d90 ,h1 ,rhosol , &
!
! k-layer contains aks (take part above)
!
if (concin(k-1)<1.0e-6_fp .or. concin(k)<1.0e-6_fp) then
ceavg = ceavg + concin(k)*(1.0_fp-dif_aks/thick(k))*thick(k)*h1
if (dif_upp <= thick(k)) then ! this must be true if k = 1 since aks < 3*deltas < h1
! k-layer contains also 3*deltas (take part below)
ceavg = ceavg + concin(k)*(dif_upp-dif_aks)*h1
exit
elseif (concin(k-1)<1.0e-6_fp .or. concin(k)<1.0e-6_fp) then
ceavg = ceavg + concin(k)*(thick(k)-dif_aks)*h1
else
rpower = log(concin(k-1)/concin(k)) / log( (h1*(1.0_fp+sig(k ))*(h1-h1*(1.0_fp+sig(k-1)))) &
& / (h1*(1.0_fp+sig(k-1))*(h1-h1*(1.0_fp+sig(k )))) )
z = ((1.0_fp+sig(k)+0.5_fp*thick(k))*h1 + aks) / 2.0_fp
ceavgtmp = concin(k) * ((h1*(1.0_fp+sig(k))*(h1-z))/(z*(h1-h1*(1.0_fp+sig(k)))))**rpower
ceavg = ceavg + ceavgtmp*(1.0_fp-dif_aks/thick(k))*thick(k)*h1
ceavg = ceavg + ceavgtmp*(thick(k)-dif_aks)*h1
endif
elseif (dif_aks<=0.0_fp .and. dif_upp>=0.0_fp) then
!
Expand Down
2 changes: 1 addition & 1 deletion src/utils_gpl/morphology/packages/morphology_kernel/src/calseddf1993.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ subroutine calseddf1993(ustarc ,ws ,h1 ,num_layers_grid ,s
if (ltur==0 .or. ltur==1 .or. difvr) then
!
! if algebraic or K-L turbulence model or difvr = .true. then
! calculate sediment mixing according to Van Rijn based on his
! calculate sediment mixing according to Van Rijn, using Coleman's
! parabolic-linear mixing distribution for current-related mixing
!
! set vertical sediment mixing values for waves and currents at water surface
Expand Down
132 changes: 72 additions & 60 deletions src/utils_gpl/morphology/packages/morphology_kernel/src/trab19.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,15 +24,12 @@
! Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!
!

!> computes sediment transport according to the transport formula of Van Thiel / Van Rijn (2008)
subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h ,tp , &
& d50 ,d15 ,d90 ,npar ,par ,dzbdt ,vicmol , &
& poros ,chezy ,dzdx ,dzdy ,sbotx ,sboty ,cesus , &
& ua ,va ,ubot ,kwtur ,ubot_from_com )
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
use precision
use mathconsts
Expand All @@ -49,17 +46,17 @@ subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h
real(fp) , intent(in) :: d50
real(fp) , intent(in) :: d90
real(fp) , intent(in) :: dzbdt !< Erosion/sedimentation velocity
real(fp) , intent(in) :: dzdx
real(fp) , intent(in) :: dzdy
real(fp) , intent(in) :: dzdx
real(fp) , intent(in) :: dzdy
real(fp) :: h
real(fp) :: hrms
real(fp) , intent(in) :: kwtur !< Breaker induced turbulence
real(fp), dimension(npar), intent(in) :: par
real(fp) , intent(in) :: poros
real(fp) , intent(in) :: rlabda
real(fp) , intent(in) :: teta
real(fp) :: tp
real(fp) , intent(in) :: ubot
real(fp) , intent(in) :: rlabda
real(fp) , intent(in) :: teta
real(fp) :: tp
real(fp) , intent(in) :: ubot
real(fp) , intent(in) :: u
real(fp) , intent(in) :: v
real(fp) , intent(in) :: vicmol
Expand All @@ -74,6 +71,7 @@ subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h
!
real(fp), parameter :: DTOL = 1e-6_fp
real(fp), parameter :: ONETHIRD = 1.0_fp/3.0_fp
real(fp), parameter :: D50_REF = 0.000225_fp !< Reference sand diameter [m]

integer :: waveform
integer :: dilatancy
Expand All @@ -90,19 +88,20 @@ subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h
real(fp) :: cmax
real(fp) :: reposeangle
real(fp) :: rheea
real(fp) :: cf
real(fp) :: utot
real(fp) :: uamag
real(fp) :: phi
real(fp) :: b2
real(fp) :: ucrw
real(fp) :: ucrc
real(fp) :: dster
real(fp) :: cf
real(fp) :: utot !< Velocity magnitude
real(fp) :: uamag
real(fp) :: phi
real(fp) :: b2
real(fp) :: ucrw
real(fp) :: ucrc
real(fp) :: dster
real(fp) :: ucr
real(fp) :: urms
real(fp) :: urms
real(fp) :: urms2
real(fp) :: ucrb, ucrs, asb, ass, term1, ceqb, ceqs
real(fp) :: cmax2h
real(fp) :: alfad50
!
!
!! executable statements -------------------------------------------------------
Expand All @@ -119,11 +118,17 @@ subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h
if ( utot < DTOL .or. h > 200.0_fp .or. h < 0.01_fp ) return
!
! Initialisations
!
ag = par(1)
delta = par(4)
facua = par(11)
facas = par(12)
facsk = par(13)
if (comparereal(facua, 0.0_fp, 1e-10_fp) == 0) then
facas = par(12)
facsk = par(13)
else
facas = facua
facsk = facua
end if
waveform = int(par(14))
sws = int(par(15))
lws = int(par(16))
Expand All @@ -134,10 +139,10 @@ subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h
smax = par(21)
reposeangle = par(22)
cmax = par(23)
alfad50 = par(24)
!
! limit input parameters to sensible values
!
facua = max(min(facua,1.0_fp),0.0_fp)
facas = max(min(facas,1.0_fp),0.0_fp)
facsk = max(min(facsk,1.0_fp),0.0_fp)
if (.not. (waveform==1 .or. waveform==2)) waveform=2 ! van Thiel default
Expand All @@ -146,11 +151,12 @@ subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h
if (.not. (dilatancy==1 .or. dilatancy==0)) dilatancy = 0 ! default off
rheea = max(min(rheea,2.0_fp),0.75_fp)
pormax = max(min(pormax,0.6_fp),poros)
if (.not. (bedslpeffini==0 .or. bedslpeffini==1 .or. bedslpeffini==2)) bedslpeffini=0
if (.not. (bedslpeffini==0 .or. bedslpeffini==1 .or. bedslpeffini==2)) bedslpeffini=0
smax = max(min(smax,3.0_fp),-1.0_fp)
if (smax<0.0_fp) smax=huge(0.0_fp)*1.0e-20_fp
reposeangle = max(min(reposeangle,45.0_fp),30.0_fp)
cmax = max(min(cmax,1.0_fp),0.0_fp)
alfad50 = max(min(alfad50,1.5_fp),0.0_fp)
!
cf = ag / chezy / chezy
!
Expand All @@ -166,48 +172,54 @@ subroutine trab19(u ,v ,hrms ,rlabda ,teta ,h
dster=(delta*ag/1e-12_fp)**onethird*d50 ! 1e-12 = nu**2
!
if(d50<=0.0005_fp) then
Ucrc=0.19_fp*d50**0.1_fp*log10(4.0_fp*h/d90) !Shields
Ucrw=0.24_fp*(delta*ag)**0.66_fp*d50**0.33_fp*tp**0.33_fp !Komar and Miller (1975)
Ucrc=0.19_fp*d50**0.1_fp*log10(4.0_fp*h/d90) ! Shields
Ucrw=0.24_fp*(delta*ag)**0.66_fp*d50**0.33_fp*tp**0.33_fp ! Komar and Miller (1975)
else if(d50<=0.002_fp) then
Ucrc=8.5_fp*d50**0.6_fp*log10(4.0_fp*h/d90) !Shields
Ucrw=0.95_fp*(delta*ag)**0.57_fp*d50**0.43_fp*tp**0.14_fp !Komar and Miller (1975)
Ucrc=8.5_fp*d50**0.6_fp*log10(4.0_fp*h/d90) ! Shields
Ucrw=0.95_fp*(delta*ag)**0.57_fp*d50**0.43_fp*tp**0.14_fp ! Komar and Miller (1975)
else if(d50>0.002_fp) then
Ucrc=1.3_fp*sqrt(delta*ag*d50)*(h/d50)**(0.5_fp*onethird) !Maynord (1978) --> also Neill (1968) where 1.3_fp = 1.4_fp
Ucrw=0.95_fp*(delta*ag)**0.57_fp*d50**0.43_fp*tp**0.14_fp !Komar and Miller (1975)
Ucrc=1.3_fp*sqrt(delta*ag*d50)*(h/d50)**(0.5_fp*onethird) ! Maynord (1978) --> also Neill (1968) where 1.3_fp = 1.4_fp
Ucrw=0.95_fp*(delta*ag)**0.57_fp*d50**0.43_fp*tp**0.14_fp ! Komar and Miller (1975)
end if
B2 = utot/max(utot+sqrt(urms2),5e-3_fp)
Ucr = B2*Ucrc + (1.0_fp-B2)*Ucrw !Van Rijn 2007 (Bed load transport paper)
Ucr = B2*Ucrc + (1.0_fp-B2)*Ucrw ! Van Rijn 2007 (Bed load transport paper)
!
call calculate_critical_velocities(dilatancy, bedslpeffini, dzbdt, ag, vicmol, d15, poros, pormax, rheea, delta, u, v, &
dzdx, dzdy, dtol, phi, ucr, ucrb, Ucrs)
!
! transport parameters
Asb=0.015_fp*h*(d50/h)**1.2_fp/(delta*ag*d50)**0.75_fp !bed load coefficent
Ass=0.012_fp*d50*dster**(-0.6_fp)/(delta*ag*d50)**1.2_fp !suspended load coeffient
!
! Van Rijn use Peak orbital flow velocity --> 0.64 corresponds to 0.4 coefficient regular waves Van Rijn (2007)
term1=utot**2+0.64_fp*sws*urms2
! reduce sediment suspensions for (inundation) overwash conditions with critical flow velocities
term1=min(term1,smax*ag/max(cf,1e-10_fp)*d50*delta)
term1=sqrt(term1)
!
ceqb = 0.0_fp !initialize ceqb
ceqs = 0.0_fp !initialize ceqs
!
if(term1>Ucrb .and. h>dtol) then
ceqb=Asb*(term1-Ucrb)**1.5_fp
end if
if(term1>Ucrs .and. h>dtol) then
ceqs=Ass*(term1-Ucrs)**2.4_fp
end if
!
cmax2h = cmax*h/2.0_fp
ceqb = min(ceqb, cmax2h) ! maximum equilibrium bed concentration
cesus = min(ceqs, cmax2h)/h ! m2/s/m*s/m = [-], and times rhosol in eqtran
ua = uamag*cos(teta*degrad)
va = uamag*sin(teta*degrad)
sbotx = (u+ua)*ceqb ! m2/s
sboty = (v+va)*ceqb
!
999 continue
!
! transport parameters
Asb=0.015_fp*h*(d50/h)**1.2_fp/(delta*ag*d50)**0.75_fp ! bed load coefficient
Ass=0.012_fp*d50*dster**(-0.6_fp)/(delta*ag*d50)**1.2_fp ! suspended load coefficient
!
! Van Rijn use Peak orbital flow velocity --> 0.64 corresponds to 0.4 coefficient regular waves Van Rijn (2007)
term1=utot**2+0.64_fp*sws*urms2
! reduce sediment suspensions for (inundation) overwash conditions with critical flow velocities
term1=min(term1,smax*ag/max(cf,1e-10_fp)*d50*delta)
term1=sqrt(term1)
!
if(term1>Ucrb .and. h>dtol) then
ceqb = Asb*(term1-Ucrb)**1.5_fp
else
ceqb = 0.0_fp
end if
if(term1>Ucrs .and. h>dtol) then
ceqs = Ass*(term1-Ucrs)**2.4_fp
else
ceqs = 0.0_fp
end if
!
if (alfad50 > 0.0_fp) then
ceqb = ceqb * (D50_REF/d50)**alfad50
ceqs = ceqs * (D50_REF/d50)**alfad50
endif
!
cmax2h = cmax*h/2.0_fp
ceqb = min(ceqb, cmax2h) ! maximum equilibrium bed concentration
cesus = min(ceqs, cmax2h)/h ! maximum equilibrium suspended concentration [m2/s/m*s/m] = [-], and times rhosol in eqtran
ua = uamag*cos(teta*degrad)
va = uamag*sin(teta*degrad)
sbotx = (u+ua)*ceqb ! [m2/s]
sboty = (v+va)*ceqb
!
999 continue
end subroutine trab19
Loading
Loading