module sandmud_module ! ! Module for the description of erosion and sedimentation processes due to sand-mud ! interaction for cohesive and uncohesive regime. Erosion and deposition flxes are ! expressed respectively in [m3/m2/s] (excluding pores) and [m/s] ! contains ! !===================================================================================! subroutine sandmud_erodep(s,par,local_wetz) ! use params use spaceparams use vegetation_module ! implicit none ! ! type(parameters) :: par type(spacepars) :: s real*8, dimension(par%ngd) :: ws ! [m/s] fall velocity real*8, dimension(par%nx+1,par%ny+1) :: taub ! [N/m2] bottom shear stress real*8, dimension(par%nx+1,par%ny+1) :: U ! [m/s] velocity magnitude at cell center integer, dimension(par%nx+1,par%ny+1) :: local_wetz ! [-] updated mask of wet cells integer :: i,j ! ! ! Initialize variables ! s%ero = 0.d0 s%depo_ex = 0.d0 s%ceqbg = 0.d0 ! ! Compute fall velocity ! TODO: find a better version for cohesive sediment ! call settling(s, par, ws) ! ! Compute velocity magnitude at each cell ! U = sqrt(s%ue**2+s%ve**2) do i=1,par%nx do j=1,par%ny if (U(i,j)>=1.5) then write(*,*) 'high U velocity...' end if end do end do ! ! Compute bottom shear stress at each cell ! TODO: add effect of turbulence due to wave breaking ? ! call tau_bed(s, par, U, taub, local_wetz) ! ! Determine erosion and deposition fluxes ! call compute_fluxes( s, par, ws, taub, U ) ! ! end subroutine sandmud_erodep !===================================================================================! subroutine compute_fluxes( s, par, ws, taub, U ) ! ! Subroutine to compute erosion and deposition fluxes ! !use precision use params use spaceparams use soil_vegetation_module ! implicit none ! ! type(parameters) :: par type(spacepars) :: s real*8, dimension(par%ngd), intent(in) :: ws ! [m/s] settling velocity real*8, dimension(par%nx+1,par%ny+1), intent(in) :: taub ! [N/m2] bottom shear stress at each cell real*8, dimension(par%nx+1,par%ny+1), intent(in) :: U ! [m/s] velocity magnitude real*8, dimension(par%nx+1,par%ny+1,par%ngd) :: sour ! [m3/m2/s] erosion flux per fraction real*8, dimension(par%nx+1,par%ny+1,par%ngd) :: sink ! [m/s] deposition lfux per fraction real*8, dimension(par%nx+1,par%ny+1) :: pbbed_mud ! [-] sum of mud sediment fractions in the upper layer real*8, dimension(par%nx+1,par%ny+1) :: Mnc ! [kg/m2/s] erosion parameter non-cohesive mixtures real*8, dimension(par%nx+1,par%ny+1) :: Mc ! [kg/m2/s] erosion parameter cohesive mixtures real*8, dimension(par%nx+1,par%ny+1) :: tau_nc ! [N/m2] critical shear stress non-cohesive mixture [N/m2] real*8, dimension(par%nx+1,par%ny+1) :: tau_c ! [N/m2] critical shear stress cohesive mixture [N/m2] real*8, dimension(par%nx+1,par%ny+1) :: hloc real*8, dimension(par%nx+1,par%ny+1) :: Tsg ! [s] reference settling time for sand fractions real*8 :: depfac real*8 :: taudep ! [N/m2] critical shear stress for deposition real*8, dimension(par%nx+1,par%ny+1) :: sour_star real*8, dimension(par%nx+1,par%ny+1,par%ngd) :: max_ero integer :: i,j,jg !real*16 :: mystart,myend ! ! hloc = max(s%hh,par%eps) Mnc = 0.0d0 Mc = 0.0d0 tau_nc = 0.0d0 tau_c = 0.0d0 sour = 0.0d0 sink = 0.0d0 sour_star = 0.0d0 s%tau_bottom = taub ! ! Determine mud fraction in the top layer ! call mud_fraction( s, par, pbbed_mud) ! ! Compute change in sediment composition (e.g. based on available fractions and sediment availability) ! if (par%ngd==1) then ! ! One sediment fraction ! else if (par%ngd==2) then ! ! Two sediment fractions ! ! Determine the erosion parameters (of each cell) based on soil local composition ! call erosion_par( s , par , pbbed_mud, Mnc , Mc , & ws , tau_nc, tau_c , hloc, Tsg, & taudep, depfac, U ) s%gMnc = Mnc s%gMc = Mc ! if (par%soilveg==1) then ! ! Compute the effect of vegetation in increasing critical shear stress following ! the approach by Van der Meer et al. (2007) and Tuan and Oumeraci (2012) ! call tau_soil_veg( s, par, tau_nc, tau_c) ! end if ! ! Start to compute fluxes ! do j=1,par%ny+1 do i=1,par%nx+1 ! ! Check sediment behavior ! if (pbbed_mud(i,j)>=par%pmcr) then ! ! Erosion/deposition cohesive mixtures ! call mud_erosion( Mc(i,j), tau_c(i,j), taub(i,j), sour_star(i,j) ) ! else ! ! Erosion/deposition uncohesive mixtures ! if (par%Mnccalc==1) then ! ! Approach based on Van Ledden (2003) ! call sand_erosion_VL( Mnc(i,j), tau_nc(i,j), taub(i,j), sour_star(i,j) ) ! else ! ! Approach based on Van Rijn (1993) ! call sand_erosion_new( Mnc(i,j), tau_nc(i,j), taub(i,j), sour_star(i,j) ) ! end if end if end do end do do jg=1,par%ngd if (par%sedtype(jg)==2) then do j=1,par%ny+1 do i=1,par%nx+1 ! ! Erosion flux mud [kg/m2/s] ! sour(i,j,jg) = pbbed_mud(i,j)*sour_star(i,j) ! ! Deposition flux mud [m/s] ! sink(i,j,jg) = mud_deposition( par, s%ccg(i,j,jg), ws(jg), taub(i,j), taudep) ! end do end do else do j=1,par%ny+1 do i=1,par%nx+1 ! ! Erosion flux sand [kg/m2/s] ! sour(i,j,jg) = (1.0d0-pbbed_mud(i,j))*sour_star(i,j) ! ! Deposition flux sand [m/s] ! if (par%Mnccalc==1) then sink(i,j,jg) = depfac*ws(jg) else !sink(i,j,jg) = s%hh(i,j)/Tsg(i,j) sink(i,j,jg) = depfac*ws(jg) end if ! end do end do end if end do else ! ! More than two sediment fractions ! end if ! ! Assign erosion/deposition fluxes (excluding pores; deposition flux represents a ! settling velocity and it needs to be multiplied by sediment concetration to obtain ! the sediment deposition flux) ! s%ero = sour/par%rhos s%depo_ex = sink ! !call allowed_erosion(s,par,max_ero) !do jg=1,par%ngd ! do j=1,par%ny+1 ! do i=1,par%nx+1 ! if (s%ero(i,j,jg)) then ! write(*,*), 'check erosion' ! end if ! end do ! end do !end do ! ! end subroutine compute_fluxes !===================================================================================! subroutine settling(s, par, ws) ! ! Subroutine to determine settling velocity based on sediment size and typology, ! based on Van Rijn (1993) ! use params use spaceparams ! implicit none ! type(parameters) :: par type(spacepars) :: s real*8, dimension(par%ngd), intent(out) :: ws ! [m/s] falling velocity real*8 :: sd ! [-] "submerged" sediment density real*8 :: nu ! [m2/s] water kinematic viscosity integer :: jg ! ! nu = 1.0040E-6 do jg = 1,par%ngd sd = (par%rhos-par%rho)/par%rho if (s%d50(jg)<65.0E-6) then ws(jg) = 5.0E-04 elseif (s%d50(jg)>65.0E-6 .and. s%d50(jg)<100.0E-6) then ws(jg) = (sd-1.0d0)*par%g*s%d50(jg)**2/(18.0d0*nu) elseif (s%d50(jg)>=100.0E-6 .and. s%d50(jg)<1000.0E-6) then ws(jg) = 10.0d0*nu/s%d50(jg)*((1.0d0+0.01d0*(sd-1.0d0)*par%g*s%d50(jg)**3.0d0/nu**2.0d0)**0.5d0-1) else ws(jg) = 1.1d0*((sd-1.0d0)*par%g*s%d50(jg))**0.5d0 end if end do ! ! end subroutine settling !===================================================================================! subroutine tau_bed(s, par, U, taub, local_wetz) ! ! Subroutine to determine the bed shear stress due to waves and currents ! based on Soulsby (1997) ! use params use paramsconst use spaceparams ! implicit none type(parameters) :: par type(spacepars) :: s real*8, dimension(par%nx+1,par%ny+1) :: U ! [m/s] near bed current velocity real*8, dimension(par%nx+1,par%ny+1), intent(out) :: taub ! [N/m2] total bed shear stress real*8, dimension(par%nx+1,par%ny+1) :: fw ! [-] wave friction parameter real*8, dimension(par%nx+1,par%ny+1) :: A ! [m] semi-orbital excursion real*8, dimension(par%nx+1,par%ny+1) :: Cd ! [-] drag coefficient real*8, dimension(par%nx+1,par%ny+1) :: tau_w ! [N/m2] wave bed shear stress real*8, dimension(par%nx+1,par%ny+1) :: tau_c ! [N/m2] current bed shear stress real*8, dimension(par%nx+1,par%ny+1) :: tau_m ! [N/m2] mean shear stress due to wave-current interaction real*8 :: phi ! [rad] angle between wave direction and current direction real*8 :: eq_d50 ! [m] equivalent grain roughness to compute the friction factor real*8, dimension(par%nx+1,par%ny+1) :: hloc real*8, dimension(par%nx+1,par%ny+1) :: urms2 ! [m2/s2] square of the mean sqared root velocity real*8 :: ML real*8 :: dcfin integer, dimension(par%nx+1,par%ny+1), intent(in) :: local_wetz ! [-] updated flag for wet cells integer :: i,j,jg ! ! hloc = max(s%hh,par%eps) phi = 0.0d0 urms2 = s%urms**2 ! do jg=1,par%ngd if (par%sedtype(jg)==1) then ! (z0) in Soulsby (1997) ! This is calculated using sand grain size eq_d50 = s%d50(jg)/12 ! eq_d50 = 5.25E-6 end if end do ! ! Calculate wave friction parameter (Soulsby, 1997) ! do j=1,par%ny+1 do i=1,par%nx+1 if (s%urms(i,j)<=0.005d0) then fw(i,j) = 0.0d0 else !fw(i,j) = 1.39d0*(s%urms(i,j)*par%Trep/(2*par%px*eq_d50))**(-0.52d0) fw(i,j) = 1.39d0*(s%urms(i,j)*par%Trep/(2*par%px*par%z0))**(-0.52d0) end if end do end do ! if (par%swave==1) then ! ! Consider near be induced turbulence due to short waves ! (perhaps in this case it is better to avoid it since urms ! is just increased due to the effct of reflection) ! do j=1,par%ny+1 do i=1,par%nx+1 ! ! Compute mixing length ! ML = sqrt(s%R(i,j)*par%Trep/(par%rho*s%c(i,j))) ML = min(ML,hloc(i,j)) ! ! Exponential decay turbulence over depth ! dcfin = exp(min(100.0d0,hloc(i,j)/max(ML,0.01d0))) dcfin = min(1.0d0,1.0d0/(dcfin-1.0d0)) s%kb(i,j) = s%kturb(i,j)*dcfin if (par%turb==TURB_BORE_AVERAGED) then s%kb(i,j) = s%kb(i,j)*par%Trep/s%Tbore(i,j) end if end do end do ! ! Increase urms^2 due to turbulence ! urms2 = s%urms**2+1.45d0*s%kb end if ! if (par%lws==0) then ! ! Switch to include long wave stirring (TODO: check if implementation is ok) ! U = (1.d0/par%cats/par%Trep*par%dt)*U ! endif ! ! Determine shear stress due to waves, current and waves-current interaction ! tau_w = 0.5d0*par%rho*fw*urms2*local_wetz ! since wave field is not solved every time step it is possible ! to have wave height (and then urms) differnt from zero at dry cells Cd = (0.4d0/(log(max(hloc,10.d0*par%z0)/par%z0)-1.0d0))**2 tau_c = par%rho*Cd*U**2 ! do j=1,par%ny+1 do i=1,par%nx+1 if (tau_w(i,j)+tau_c(i,j)==0.0d0 .or. hloc(i,j)<=par%eps) then tau_m(i,j) = 0.0d0 else tau_m(i,j) = tau_c(i,j)*(1.0d0+1.2d0*(tau_w(i,j)/(tau_w(i,j)+tau_c(i,j)))**3.2d0) end if end do end do ! taub = sqrt((tau_m+tau_w*cos(phi))**2+(tau_w*sin(phi))**2) !do j=1,par%ny+1 ! do i=1,par%nx+1 ! if (taub(i,j)>=15.0d0) then ! write(*,*), 'high tau...' ! end if ! end do !end do ! ! end subroutine tau_bed !===================================================================================! subroutine erosion_par( s , par , pm , Mnc , Mc , & ws , tau_nc, tau_c, hloc, Tsg , & taudep, depfac, U ) ! ! Subroutine to determine the erosion parameters for sand-mud mixtures, at each cell, ! based on the relative fractions of sediment ! use params use spaceparams ! implicit none ! type(parameters) :: par type(spacepars) :: s real*8, dimension(par%nx+1,par%ny+1), intent(in) :: pm ! [-] relative amount of mud sediment fractions in the cell real*8, dimension(par%ngd), intent(in) :: ws ! [m/s] settiling velocity per fraction real*8, dimension(par%nx+1,par%ny+1), intent(in) :: hloc real*8, dimension(par%nx+1,par%ny+1), intent(in) :: U real*8, dimension(par%nx+1,par%ny+1), intent(out) :: Mnc ! [kg/m2/s] erosion parameter for uncohesive mixtures real*8, dimension(par%nx+1,par%ny+1), intent(out) :: Mc ! [kg/m2/s] erosion parameter for cohesive mixtures real*8, dimension(par%nx+1,par%ny+1) :: Ms ! [kg/m2/s] erosion parameter for pure sand Ms = Mnc(pm=0) real*8, dimension(par%nx+1,par%ny+1) :: Mcr ! [kg/m2/s] erosion parameter at pm = pmcr real*8, dimension(par%nx+1,par%ny+1), intent(out) :: tau_nc ! [N/m2] critical shear stress for uncohesive mixtures real*8, dimension(par%nx+1,par%ny+1), intent(out) :: tau_c ! [N/m2] critical shear stress for cohesive mixtures real*8, dimension(par%nx+1,par%ny+1), intent(out) :: Tsg real*8, intent(out) :: taudep real*8, intent(out) :: depfac real*8, dimension(par%nx+1,par%ny+1) :: expon real*8 :: nu ! [m2/s] water kinematic viscosity real*8, dimension(par%ngd) :: Dstar ! [-] dimensionless grain size real*8 :: thetacr ! [-] dimensionless critical SHields number real*8 :: sr ! [-] relative sediment density real*8 :: alfa ! [-] parameter for entrainment uncohesive mixtures real*8 :: tau_cr_sand real*8 :: Fs ! [-] Shape factor to account for sediment distribution along water column Van Rijn (1993), range 0.1-0.03 real*8 :: ac_vr ! [m] Van Rijn reference height real*8, dimension(par%nx+1,par%ny+1) :: explore integer :: jg,i,j ! ! ! Initialize variables ! Mcr = 0.0d0 Mnc = 0.0d0 Ms = 0.0d0 Mc = 0.0d0 tau_nc = 0.0d0 tau_c = 0.0d0 Tsg = 0.0d0 Dstar = 0.0d0 thetacr = 0.0d0 taudep = par%taudep depfac = 1.0d0 nu = 1.0040E-6 sr = par%rhos/par%rho ! ! Determination of Mnc; the value of Mnc has to be multiplied by the ! sediment (dry?) density to recover the right dimension [kg/m2/s] ! do jg=1,par%ngd if (par%sedtype(jg)==1) then ! ! Detrmination of Dstar and erosion parameter for uncohesive sediment ! Dstar(jg) = s%d50(jg)*((sr-1.0d0)*par%g/nu**2)**(1.0d0/3.0d0) ! if (par%Mnccalc==1) then ! ! Van Ledden (2003) ! alfa = 0.075 ! Average value for alphab1 in Van Rijn (1993) bedload formulation Ms = par%rhos*alfa*sqrt((sr-1.0d0)*par%g*s%d50(jg))/(Dstar(jg)**0.9d0) ! TODO: REVIEW THE FORMULATION, TOO HIGH VALUES !Mnc = 0.1d0 else ! ! New formulation using Van Rijn (1993) ! Fs = 0.1d0 ! Value from Van Rijn (1993) (check) !Fs = 0.2d0 ac_vr = 2.5*s%d50(jg) ! Van Rijn reference height chosen equal to the Nikuradse roughness (Soulsby, 1997) ! ! Approach dependent on water depth ! !Tsg = par%tsfac*hloc/ws(jg) !Tsg = max(Tsg,par%Tsmin) !Tsg = max(0.1*hloc/ws(jg),0.2) !Ms = par%rhos*1.5*Fs*par%d50(jg)/(Dstar(jg)**0.3*Tsg) ! ! Approach not dependent on water depth ! Ms = par%rhos*ws(jg)*Fs*0.015*s%d50(jg)/(ac_vr*Dstar(jg)**0.3) end if ! ! Determination of the critical shear stress from Shields (1932) ! thetacr = det_shields(Dstar(jg)) tau_cr_sand = (par%rhos-par%rho)*par%g*s%d50(jg)*thetacr ! tau_cr_sand = par%tausand ! Dano: this looks fishy! ! if (tau_cr_sand<=0.0d0) then ! tau_cr_sand = par%tausand ! end if end if end do ! !if (Mnc>0.20d0 .or. Mnc<0.005) then ! ! ! ! Value of Mnc is set as suggested by Van Ledden (2003) ! ! ! Mnc = 0.05d0 ! expressed in [kg/m2/s] !end if ! ! Determination of critical erosion parameter ! Mcr = par%Mmud*par%betacr+(1.0d0-par%betacr)*Ms ! ! Determination of erosion parameter for non-cohesive mixtures ! explore = min(pm,par%pmcr) Mnc = Ms+(Mcr-Ms)/par%pmcr*explore ! ! Determination of erosion parameter for cohesive mixtures ! expon = min((1.0d0-pm)/(1.0d0-par%pmcr),1.0d0) Mc = par%Mmud*(Mcr/par%Mmud)**expon ! ! Determination of critical shear stresses ! tau_nc = tau_cr_sand*(1.0d0+pm)**par%betam tau_c = (tau_cr_sand*(1.0d0+par%pmcr)**par%betam-par%taumud)/(1.0d0-par%pmcr)*(1.0d0-pm)+par%taumud ! ! end subroutine erosion_par !===================================================================================! subroutine mud_fraction( s, par, pbbed_mud) ! ! Subroutine to get mud fraction in a cell ! use spaceparams use params ! implicit none ! ! type(spacepars) :: s type(parameters) :: par real*8, dimension(par%nx+1,par%ny+1) :: pbbed_mud integer :: i,j,jg ! ! pbbed_mud = 0.0d0 do jg=1,par%ngd if (par%sedtype(jg)==2) then do j=1,par%ny+1 do i=1,par%nx+1 pbbed_mud(i,j) = pbbed_mud(i,j)+s%pbbed(i,j,1,jg) end do end do !if (pbbed_mud(84,1)0.0d0) then taum = max(0.0d0,taub/tau_cr-1.0d0) end if ! ! Erosion/deposition flux in cohesive regime ! sour_star = Mc*taum ! ! end subroutine mud_erosion !===================================================================================! function mud_deposition( par, c_mud, ws, taub, taudep) ! ! Subroutine to determine mud deposition based on critical ! shear stress for deposition ! use params ! implicit none ! ! type(parameters) :: par real*8, intent(in) :: c_mud real*8, intent(in) :: ws real*8, intent(in) :: taub real*8, intent(in) :: taudep real*8 :: mud_deposition real*8 :: sink_star real*8 :: k_mud_dep ! [] coefficient for mud settling velocity due to flocculation real*8 :: ws_c ! [m/s] concentration dependent settling velocity real*8 :: c_mass ! [g/l] massive mud conentration ! ! k_mud_dep = 0.001d0 ws_c = 0.d0 sink_star = 0.0d0 if (taudep>0.0d0) then sink_star = max(0.0d0,1-taub/taudep) end if ! c_mass = c_mud*par%rhos*par%morfac !ws_c = max(c_mass**(4.d0/3.d0)*k_mud_dep, ws) ws_c = max(c_mass*k_mud_dep, ws) mud_deposition = ws_c*sink_star ! ! end function mud_deposition !===================================================================================! subroutine sand_erosion_VL( Mnc, tau_cr, taub, sour_star ) ! ! Van Ledden (2003) approach ! implicit none ! ! real*8, intent(in) :: Mnc real*8, intent(in) :: tau_cr real*8, intent(in) :: taub real*8, intent(out) :: sour_star real*8 :: taum ! ! sour_star = 0.0d0 taum = 0.0d0 if (tau_cr>0.0d0) then taum = max(0.0d0,taub/tau_cr-1.0d0) end if ! ! Erosion flux for sand non-cohesive regime [kg/m2/s] ! sour_star = Mnc*taum ! ! end subroutine sand_erosion_VL !===================================================================================! subroutine sand_erosion_new( Mnc, tau_cr, taub, sour_star ) ! ! Subroutine to determine erosion flux for the sand fractions ! in uncohesive regime based on Van Rijn formulation (Van Rijn, 1993). ! use params use spaceparams ! implicit none ! ! real*8, intent(in) :: Mnc real*8, intent(in) :: tau_cr real*8, intent(in) :: taub real*8, intent(out) :: sour_star real*8 :: taum ! ! sour_star = 0.0d0 taum = 0.0d0 if (tau_cr>0.0d0) then taum = max(0.0d0,taub/tau_cr-1.0d0) end if ! ! Erosion flux for sand non-cohesive regime [kg/m2/s] ! sour_star = Mnc*taum !**1.5d0 ! ! end subroutine sand_erosion_new !===================================================================================! function det_shields(Dstar) ! ! SUbroutine to determine Shields number for a ! specific sand grain size based on Van Rijn (1993) ! implicit none ! ! real*8, intent(in) :: dstar real*8 :: det_shields ! ! if (Dstar<1.0d0) then det_shields = 0.0d0 else if (Dstar<=4.0d0) then det_shields = 0.24d0/Dstar else if (Dstar<=10.0d0) then det_shields = 0.140/Dstar**0.64d0 else if (Dstar<=20.0d0) then det_shields = 0.04d0/Dstar**0.1d0 else if (Dstar<=150.0d0) then det_shields = 0.013d0*Dstar**0.29d0 else det_shields = 0.055d0 endif ! ! end function det_shields !===================================================================================!! subroutine allowed_erosion(s,par,max_ero) ! ! ! use params use spaceparams ! implicit none ! ! type(parameters) :: par type(spacepars) :: s real*8, dimension(par%nx+1,par%ny+1,par%ngd), intent(out) :: max_ero integer :: i,j,jg ! ! max_ero = 0.0d0 do jg=1,par%ngd do j=1,par%ny+1 do i=1,par%nx+1 max_ero(i,j,jg) = s%pbbed(i,j,1,jg)*s%dzbed(i,j,1)*(1.0d0-par%por)/par%dt/par%morfac end do end do end do ! ! end subroutine allowed_erosion !===================================================================================!! ! DIFFERENT WAY TO COMPUTE FLUXES !call cpu_time(mystart) !do jg=1,par%ngd ! if (par%sedtype(jg)==2) then ! do j=1,par%ny+1 ! do i=1,par%nx+1 ! ! ! ! Check sediment behavior ! ! ! if (pbbed_mud(i,j)>=par%pmcr) then ! ! ! ! Erosion/deposition cohesive mixtures ! ! ! call mud_erosion( Mc(i,j), tau_c(i,j), taub(i,j), sour_star ) ! ! ! else ! ! ! ! Erosion/deposition uncohesive mixtures ! ! ! if (par%Mnccalc==1) then ! ! ! ! Approach based on Van Ledden (2003) ! ! ! call sand_erosion_VL( Mnc(i,j), tau_nc(i,j), taub(i,j), sour_star ) ! ! ! else ! ! ! ! Approach based on Soulsby-Van Rijn ! ! ! call sand_erosion_SBVR( s , par , sour_star, U(i,j), pbbed_mud(i,j), & ! Mnc(i,j), Tsg(i,j), hloc(i,j), i , j ) ! ! ! end if ! end if ! ! ! ! Erosion flux mud [kg/m2/s] ! ! ! sour(i,j,jg) = pbbed_mud(i,j)*sour_star ! ! ! ! Deposition flux mud [m/s] ! ! ! sink(i,j,jg) = mud_deposition( ws(jg), taub(i,j), taudep) ! ! ! end do ! end do ! else ! do j=1,par%ny+1 ! do i=1,par%nx+1 ! ! ! ! Check sediment behavior ! ! ! if (pbbed_mud(i,j)>=par%pmcr) then ! ! ! ! Erosion/deposition cohesive mixtures ! ! ! call mud_erosion( Mc(i,j), tau_c(i,j), taub(i,j), sour_star ) ! ! ! else ! ! ! ! Erosion/deposition uncohesive mixtures ! ! ! if (par%Mnccalc==1) then ! ! ! ! Approach based on Van Ledden (2003) ! ! ! call sand_erosion_VL( Mnc(i,j), tau_nc(i,j), taub(i,j), sour_star ) ! ! ! else ! ! ! ! Approach based on Soulsby-Van Rijn ! ! ! call sand_erosion_SBVR( s , par , sour_star, U(i,j), pbbed_mud(i,j), & ! Mnc(i,j), Tsg(i,j), hloc(i,j), i , j ) ! ! ! end if ! end if ! ! ! ! Erosion flux sand [kg/m2/s] ! ! ! sour(i,j,jg) = (1-pbbed_mud(i,j))*sour_star ! ! ! ! Deposition flux sand [m/s] ! ! ! if (par%Mnccalc==1) then ! sink(i,j,jg) = depfac*ws(jg) ! else ! sink(i,j,jg) = s%hh(i,j)/Tsg(i,j) ! end if ! ! ! end do ! end do ! end if !end do !call cpu_time(myend) !write (*,*) 'time needed', myend-mystart end module sandmud_module