diff --git a/src/utils_gpl/morphology/packages/morphology_io/src/rdtrafrm.f90 b/src/utils_gpl/morphology/packages/morphology_io/src/rdtrafrm.f90 index 6394436193c..1424c247022 100644 --- a/src/utils_gpl/morphology/packages/morphology_io/src/rdtrafrm.f90 +++ b/src/utils_gpl/morphology/packages/morphology_io/src/rdtrafrm.f90 @@ -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' @@ -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' @@ -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 diff --git a/src/utils_gpl/morphology/packages/morphology_kernel/src/bedbc2004.f90 b/src/utils_gpl/morphology/packages/morphology_kernel/src/bedbc2004.f90 index ec191c7dee5..5c6d41e6204 100644 --- a/src/utils_gpl/morphology/packages/morphology_kernel/src/bedbc2004.f90 +++ b/src/utils_gpl/morphology/packages/morphology_kernel/src/bedbc2004.f90 @@ -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 diff --git a/src/utils_gpl/morphology/packages/morphology_kernel/src/bedtr2004.f90 b/src/utils_gpl/morphology/packages/morphology_kernel/src/bedtr2004.f90 index 807387726bb..b5fe9e23128 100644 --- a/src/utils_gpl/morphology/packages/morphology_kernel/src/bedtr2004.f90 +++ b/src/utils_gpl/morphology/packages/morphology_kernel/src/bedtr2004.f90 @@ -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 ! diff --git a/src/utils_gpl/morphology/packages/morphology_kernel/src/calseddf1993.f90 b/src/utils_gpl/morphology/packages/morphology_kernel/src/calseddf1993.f90 index a460071d44a..5ad0fadb580 100644 --- a/src/utils_gpl/morphology/packages/morphology_kernel/src/calseddf1993.f90 +++ b/src/utils_gpl/morphology/packages/morphology_kernel/src/calseddf1993.f90 @@ -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 diff --git a/src/utils_gpl/morphology/packages/morphology_kernel/src/trab19.f90 b/src/utils_gpl/morphology/packages/morphology_kernel/src/trab19.f90 index bdb57be4a0c..50fccf9f247 100644 --- a/src/utils_gpl/morphology/packages/morphology_kernel/src/trab19.f90 +++ b/src/utils_gpl/morphology/packages/morphology_kernel/src/trab19.f90 @@ -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 @@ -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 @@ -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 @@ -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 ------------------------------------------------------- @@ -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)) @@ -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 @@ -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 ! @@ -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 diff --git a/src/utils_gpl/morphology/packages/morphology_kernel/src/trab20.f90 b/src/utils_gpl/morphology/packages/morphology_kernel/src/trab20.f90 index 45f777587e0..de7d7edc896 100644 --- a/src/utils_gpl/morphology/packages/morphology_kernel/src/trab20.f90 +++ b/src/utils_gpl/morphology/packages/morphology_kernel/src/trab20.f90 @@ -24,8 +24,6 @@ ! Stichting Deltares. All rights reserved. ! !------------------------------------------------------------------------------- -! -! !> computes sediment transport according to the transport formula of Soulsby / Van Rijn, XBeach flavour subroutine trab20(u ,v ,hrms ,rlabda ,teta ,h ,tp , & @@ -48,17 +46,17 @@ subroutine trab20(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) :: 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 @@ -73,6 +71,7 @@ subroutine trab20(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 @@ -91,26 +90,29 @@ subroutine trab20(u ,v ,hrms ,rlabda ,teta ,h real(fp) :: rheea real(fp) :: cf real(fp) :: utot !< Velocity magnitude - real(fp) :: uamag - real(fp) :: phi - real(fp) :: dster - real(fp) :: ucr - real(fp) :: urms + real(fp) :: uamag + real(fp) :: phi + real(fp) :: dster + real(fp) :: ucr + real(fp) :: urms real(fp) :: urms2 real(fp) :: ucrb, ucrs, asb, ass, term1, ceqb, ceqs real(fp) :: z0 real(fp) :: cd real(fp) :: cmax2h + real(fp) :: alfad50 ! ! !! executable statements ------------------------------------------------------- ! + ! + ! Initialize Transports to zero + ! sbotx = 0.0_fp sboty = 0.0_fp - cesus = 0.0_fp ua = 0.0_fp va = 0.0_fp - + cesus = 0.0_fp utot = sqrt(u**2 + v**2) if ( utot < DTOL .or. h > 200.0_fp .or. h < 0.01_fp ) return ! @@ -119,8 +121,13 @@ subroutine trab20(u ,v ,hrms ,rlabda ,teta ,h 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)) @@ -132,10 +139,10 @@ subroutine trab20(u ,v ,hrms ,rlabda ,teta ,h reposeangle = par(22) cmax = par(23) z0 = par(24) + alfad50 = par(25) ! ! 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 @@ -144,60 +151,71 @@ subroutine trab20(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) z0 = max(min(z0,0.05_fp),0.0001_fp) + alfad50 = max(min(alfad50,1.5_fp),0.0_fp) ! cf = ag / chezy / chezy ! call calculate_urms(hrms, tp, h, ag, ubot_from_com, ubot, kwtur, urms, urms2) ! + ! velocity asymmetry + ! call calculate_velocity_asymmetry(waveform, facas, facsk, sws, h, hrms, rlabda, ag, tp, urms, uamag) ! + ! Velocity magnitude + ! phi = reposeangle*degrad ! Angle of internal friction dster=(delta*ag/1e-12_fp)**onethird*d50 ! 1e-12 = nu**2 ! if(d50<=0.0005_fp) then Ucr=0.19_fp*d50**0.1_fp*log10(4.0_fp*h/d90) - else + else Ucr=8.5_fp*d50**0.6_fp*log10(4.0_fp*h/d90) - end if + end if ! call calculate_critical_velocities(dilatancy, bedslpeffini, dzbdt, ag, vicmol, d15, poros, pormax, rheea, delta, u, v, & dzdx, dzdy, dtol, phi, ucr, ucrb, Ucrs) - ! - ! drag coefficient - Cd=(0.40_fp/(log(max(h,10.0_fp*z0)/z0)-1.0_fp))**2 - ! - ! transport parameters - Asb=0.005_fp*h*(d50/h/(delta*ag*d50))**1.2_fp ! bed load coefficent - Ass=0.012_fp*d50*dster**(-0.6_fp)/(delta*ag*d50)**1.2_fp ! suspended load coeffient - ! - term1=utot**2+0.018_fp/Cd*sws*urms2 - ! - term1=min(term1,smax*ag/cf*d50*delta) - term1=sqrt(term1) - ! - ceqb = 0.0_fp - ceqs = 0.0_fp - ! - if(term1>Ucrb .and. h>dtol) then - ceqb=Asb*(term1-Ucrb)**2.4_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 ! maximum equilibrium suspended concentration - ua = uamag*cos(teta*degrad) - va = uamag*sin(teta*degrad) - sbotx = (u+ua)*ceqb - sboty = (v+va)*ceqb - ! - 999 continue + ! + ! drag coefficient + Cd=(0.40_fp/(log(max(h,10.0_fp*z0)/z0)-1.0_fp))**2 + ! + ! transport parameters + Asb=0.005_fp*h*(d50/h/(delta*ag*d50))**1.2_fp ! bed load coefficient + Ass=0.012_fp*d50*dster**(-0.6_fp)/(delta*ag*d50)**1.2_fp ! suspended load coefficient + ! + term1=utot**2+0.018_fp/Cd*sws*urms2 + ! + term1=min(term1,smax*ag/cf*d50*delta) + term1=sqrt(term1) + ! + if(term1>Ucrb .and. h>dtol) then + ceqb = Asb*(term1-Ucrb)**2.4_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 trab20