subroutine sandmud_sb_vr(roller   , phase_vel , ws      , ue    , ve         , &
                       & bed_turb , turb      , bore    , urms  , fric_coeff , &
                       & eqcs     , eqcb      , rhosoil , depth , myd50      , &
                       & myd90)
                       
    use params
    
    implicit none
    type(params)         :: par
    real*8 , intent(in)  :: roller              
    real*8 , intent(in)  :: phase_vel 
    real*8 , intent(in)  :: ws
    real*8 , intent(in)  :: ue 
    real*8 , intent(in)  :: ve
    real*8 , intent(in)  :: bed_turb              
    real*8 , intent(in)  :: turb 
    real*8 , intent(in)  :: bore
    real*8 , intent(in)  :: urms
    real*8 , intent(in)  :: fric_coeff
    real*8 , intent(in)  :: rhosoil
    real*8 , intent(in)  :: depth
    real*8 , intent(in)  :: myd50
    real*8 , intent(in)  :: myd90
    real*8 , intent(out) :: eqcs
    real*8 , intent(out) :: eqcb
    ! local
    real*8               :: ML
    real*8               :: hloc
    real*8               :: dcfin
    real*8               :: dcf
    real*8               :: wmg
    real*8               :: urms2
    real*8               :: Ts
    real*8               :: Ucr
    !
    hloc = min(depth,0.01)
    vmg  = 0.d0
    !
    ! Wave breaking induced turbulence due to short waves
    ! Compute mixing length
    !
    ML    = sqrt(2*roller*par%Trep/(rhosoil*phase_vel))
    ML    = min(ML,hloc)
    dcfin = exp(min(100.d0,hloc/max(ML,0.01d0)))
    dcf   = min(1.d0,1.d0/(dcfin-1.d0))
    !
    ! Turbulence
    !
    bed_turb = turb*dcf
    if (trim(par%turb)=='bore_averaged') then
       bed_turb = bed_turb*par%Trep/bore
    end if
    !
    ! Including long wave stirring
    !
    if (par%lws==1) then
       vmg = sqrt(ue**2+ve**2)
    elseif (par%lws==0) then
       vmg = (1.d0-1.d0/par%cats/par%Trep*par%dt)*vmg + (1.d0/par%cats/par%Trep*par%dt)*sqrt(ue**2+ve**2) 
    end if
    !
    urms2 = urms**2+1.45d0*bed_turb
    !
    Ts = par%tsfac*hloc/ws
    Ts = max(Ts,par%Tsmin)
    !
    ! Determination of critical velocity    
    !
    if (myd50<=0.0005d0) then
       Ucr = 0.19d0*myd50**0.1d0*log10(4.d0*hloc/myd90)
    else
    
    
    
    
end subroutine sandmud_sb_vr