subroutine transus(s,par)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University !
! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski,            !
! Jaap van Thiel de Vries, Robert McCall                                  !       
!                                                                         !
! d.roelvink@unesco-ihe.org                                               !
! UNESCO-IHE Institute for Water Education                                !
! P.O. Box 3015                                                           !
! 2601 DA Delft                                                           !
! The Netherlands                                                         !
!                                                                         !
! This library is free software; you can redistribute it and/or           !
! modify it under the terms of the GNU Lesser General Public              !
! License as published by the Free Software Foundation; either            !
! version 2.1 of the License, or (at your option) any later version.      !
!                                                                         !
! This library 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        !
! Lesser General Public License for more details.                         !
!                                                                         !
! You should have received a copy of the GNU Lesser General Public        !
! License along with this library; if not, write to the Free Software     !
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307     !
! USA                                                                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use params
use spaceparams

IMPLICIT NONE

type(spacepars),target              :: s
type(parameters)                    :: par

integer                             :: i
integer                             :: j,jg
real*8,dimension(:,:),allocatable,save   :: vmag2
real*8 , dimension(s%nx+1,s%ny+1)   :: dzbx,dzby
real*8,dimension(:,:),allocatable,save   :: cc,cu,cv,Su,Sv,Dc,termh

include 's.ind'
include 's.inp'

if (.not. allocated(vmag2)) then
   allocate(vmag2 (s%nx+1,s%ny+1))
   allocate(cu (s%nx+1,s%ny+1))
   allocate(cv (s%nx+1,s%ny+1))
   allocate(cc (s%nx+1,s%ny+1))
   allocate(Su (s%nx+1,s%ny+1))
   allocate(Sv (s%nx+1,s%ny+1))
   allocate(Dc (s%nx+1,s%ny+1))
   allocate(termh (s%nx+1,s%ny+1))
endif
! use eulerian velocities
vmag2    = ue**2+ve**2;
dcdx     = 0.0d0
dcdy     = 0.0d0
! calculate equilibrium concentration
! soulsby van Rijn
call sb_vr(s,par)

! compute diffusion coefficient

if (par%nuhfac==1) then
   termh = hh/max(H,.01)
   ! termh = max(hh(i,j)/2.d0,0.01d0)
   termh = min(termh,10.);
   ! Dc = 5*(DR/par%rho)**(1.d0/3.d0)*hh/(exp(termh)-1.d0)
   ! Dc = 5*par%nuhfac*hh*(DR/par%rho)**(1.d0/3.d0)
   Dc = par%dico
else
   Dc = par%dico
end if

do jg = 1,par%ngd
   cc = ccg(:,:,jg)
   !
   ! X-direction
   !
   do j=1,ny+1
      do i=1,nx
         if(ueu(i,j)>0.d0) then
            cu(i,j)=cc(i,j)
         elseif(ueu(i,j)<0.d0) then
            cu(i,j)=cc(i+1,j)
         else
            cu(i,j)=0.5d0*(cc(i,j)+cc(i+1,j))
         endif
         dcdx(i,j)=(cc(i+1,j)-cc(i,j))/(xz(i+1)-xz(i)) !Jaap
      enddo 
   enddo
   cu(nx+1,:) = cc(nx+1,:) !Robert

    !
    ! Bed slope terms
    !
    dzbx=0.d0
    do j=1,ny+1
        do i=1,nx
            dzbx(i,j)=(zb(i+1,j)-zb(i,j))/(xz(i+1)-xz(i))
        enddo
    enddo
    Su=(cu*ueu*hu-Dc*hu*dcdx-par%facsl*cu*vmagu*hu*dzbx)*wetu   !

    !Su=(cu*ueu*hu-Dc*hu*dcdx)*wetu !
!   if (par%t>=85)then
!      Su=Su
!   end if
   !
   ! Y-direction
   !
   do j=1,ny
      do i=1,nx+1
         if(vev(i,j)>0) then
            cv(i,j)=cc(i,j)
         else if(vev(i,j)<0) then
            cv(i,j)=cc(i,j+1)
         else
         cv(i,j)=0.5d0*(cv(i,j)+cv(i,j+1))
         end if
         dcdy(i,j)=(cc(i,j+1)-cc(i,j))/(yz(j+1)-yz(j)) !Jaap
      end do 
    end do
    cv(:,ny+1) = cc(:,ny+1) !Robert


    !
    ! Bed slope terms
    !
    dzby=0.d0
    do j=1,ny
        do i=1,nx+1
            dzby(i,j)=(zb(i,j+1)-zb(i,j))/(yz(j+1)-yz(j))
        enddo
    enddo
    Sv=(cv*vev*hv-Dc*hv*dcdy-par%facsl*cv*vmagv*hv*dzby)*wetv
    !Sv=(cv*vev*hv-Dc*hv*dcdy)*wetv
    
    do j=2,ny+1
       do i=2,nx+1
          cc(i,j) = hold(i,j)*cc(i,j)-par%dt*((Su(i,j)-Su(i-1,j))/(xu(i)-xu(i-1))+&
                                              (Sv(i,j)-Sv(i,j-1))/(yv(j)-yv(j-1))-&
                                               hold(i,j)*(ceqg(i,j,jg)*graindistr(i,j,1,jg)-cc(i,j))/Tsg(i,j,jg)) ! Jaap: use Tsg instead of Ts
                                             ! hold(i,j)*(ceqg(i,j,jg)*graindistr(i,j,1,jg)-cc(i,j))/Ts(i,j)) 

          cc(i,j)=max(cc(i,j),0.0d0)
       enddo
    enddo

    do j=1,ny+1
       do i=1,nx+1
          if(hh(i,j)>=par%hmin) then !Jaap
             cc(i,j)=cc(i,j)/hh(i,j)    !Jaap
          else
             cc(i,j)=0.d0
          end if
       end do
    end do
    cc(1,:)=cc(2,:)
    cc(:,1)=cc(:,2)
    cc(nx+1,:)=cc(nx+1-1,:)
    cc(:,ny+1)=cc(:,ny+1-1)
    !Jaap
    cc=cc*wetz
    !
    ccg(:,:,jg) = cc
    Svg(:,:,jg) = Sv
    Sug(:,:,jg) = Su
end do

vmag=sqrt(max(vmag2,par%umin))

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine bed_update(s,par)
use params
use spaceparams

IMPLICIT NONE

type(spacepars),target              :: s
type(parameters)                    :: par

integer                                   :: i,ii
integer                                   :: j,jg,jd
real*8                                    :: dzb,dzmax
real*8 , dimension(s%nx+1,s%ny+1)         :: dzbx,dzby,hloc,Su,Sv 
real*8 , dimension(s%nx+1,s%ny+1)         :: dzbtot,fact0,fact1,fact2,fact3,fact4,fact5,dry
real*8 , dimension(s%nx+1,s%ny+1,par%ngd) :: dzbdtg
real*8 , dimension(s%nx+1,s%ny+1,par%nd)  :: graindistrm

include 's.ind'
include 's.inp'

dzbdtg = 0.d0
dzbtot = 0.d0

! Update bed level using continuity eq.
dzb=0.d0 !!!Ap
dzbdt=0.0d0

if (par%t>=par%morstart) then

do jg = 1,par%ngd
   ! Update bed level using continuity eq.
   Su = Sug(:,:,jg)
   Sv = Svg(:,:,jg)

   ! Ap and Jaap: Bottom update is morfac dependent....
   do ii=1,par%morfac   
      do j=2,ny !Jaap nx instead of ny+1
         do i=2,nx  !Jaap nx instead of nx+1
            dzbdt(i,j)=1/(1-par%por)*        &
                (-(Su(i,j)-Su(i-1,j))/(xu(i)-xu(i-1)) &
                 -(Sv(i,j)-Sv(i,j-1))/(yv(j)-yv(j-1)))
            zb(i,j)=zb(i,j)+dzbdt(i,j)*par%dt
            sedero(i,j)=sedero(i,j)+dzbdt(i,j)*par%dt           ! total bed level changes for all sediment clases over all time steps so far
            dzbdtg(i,j,jg) = dzbdtg(i,j,jg) + dzbdt(i,j)*par%dt ! total bed level change for each sediment class at this time step
            dzbtot(i,j) = dzbtot(i,j) +  dzbdt(i,j)*par%dt      ! total bed level change for all sediment classes at this time step
          enddo
       enddo
    enddo
enddo

fact0 = 0.d0
fact1 = 0.d0
fact2 = 0.d0
fact3 = 0.d0
fact4 = 0.d0
fact5 = 0.d0

! update sediment fractions in sediment layers
do j=2,ny 
   do i=2,nx  
      if (dzbtot(i,j)<0.d0) then
         fact0(i,j) = 0.d0      
         fact1(i,j) = 1.d0
         fact2(i,j) = -dzbtot(i,j)/par%dzg
         fact3(i,j) = 0.d0
         fact4(i,j) = (par%dzg+dzbtot(i,j))/par%dzg
         fact5(i,j) = -dzbtot(i,j)/par%dzg
      else
         fact0(i,j) = 1.d0  ! only in case of vegetated layer 
         fact1(i,j) = (par%dzg-dzbtot(i,j))/par%dzg
         fact2(i,j) = 0.d0
         fact3(i,j) = dzbtot(i,j)/par%dzg
         fact4(i,j) = (par%dzg-dzbtot(i,j))/par%dzg
         fact5(i,j) = 0.d0
      endif
   enddo
enddo

do jd = 1,par%nd-1
   do jg = 1,par%ngd
      if (jd == 1) then
         if (jg == 1) then
            graindistr(:,:,jd,jg) = dzbdtg(:,:,jg)/par%dzg + fact1*graindistr(:,:,jd,jg) + fact2*graindistr(:,:,jd+1,jg) + max(dzbdtg(:,:,jg+1)/par%dzg,0.)
         elseif (jg==2) then
            !graindistr(:,:,jd,jg) = dzbdtg(:,:,jg)/par%dzg + fact1*graindistr(:,:,jd,jg) + fact2*graindistr(:,:,jd+1,jg)
            graindistr(:,:,jd,jg) = min(dzbdtg(:,:,jg)/par%dzg,0.) + fact1*graindistr(:,:,jd,jg) + fact2*graindistr(:,:,jd+1,jg)
         else
            graindistr(:,:,jd,jg) = dzbdtg(:,:,jg)/par%dzg + fact1*graindistr(:,:,jd,jg) + fact2*graindistr(:,:,jd+1,jg) 
         endif
      else
         graindistr(:,:,jd,jg) = fact3*graindistr(:,:,jd-1,jg) + fact4*graindistr(:,:,jd,jg) + fact5*graindistr(:,:,jd+1,jg)
      endif
   enddo
enddo

! ensure total fractions equal 1
graindistrm = 0.d0
do jd = 1,par%nd
    do jg = 1,par%ngd
       graindistrm(:,:,jd) = graindistrm(:,:,jd) + graindistr(:,:,jd,jg)
    enddo
    do jg = 1,par%ngd
       graindistr(:,:,jd,jg) = graindistr(:,:,jd,jg)/max(graindistrm(:,:,jd),.001)
    enddo
enddo
!
! bed boundary conditions
! 
! Fix bed at back boundary, but allow sediment to pass though
!zb(nx+1,:) = zb(nx+1,:)+dzbdt(nx,:)*par%dt  !Ap
!zb(nx+1,:) = zb(nx+1,:)+dzbdt(nx,:)*par%dt  !Ap
!sedero(nx+1,:) = sedero(nx,:) !Ap
zb(:,1) = zb(:,2)
sedero(:,1) = sedero(:,2)
zb(:,ny+1) = zb(:,ny)
sedero(:,ny+1) = sedero(:,ny)
!
! Avalanching
!
do ii=1,par%morfac
   dzbx=0*zb
   dzby=0*zb
   do j=1,ny+1
      do i=1,nx
         dzbx(i,j)=(zb(i+1,j)-zb(i,j))/(xz(i+1)-xz(i))
      enddo
   enddo
   do j=1,ny
      do i=1,nx+1
         dzby(i,j)=(zb(i,j+1)-zb(i,j))/(yz(j+1)-yz(j))
      enddo
   enddo
   !
   dry = abs(1-wetu)
   do i=2,nx-1
      do j=1,ny+1
         if(max(hh(i,j),hh(i+1,j))>par%hswitch+par%eps) then
            dzmax=par%wetslp;
         else
            dzmax=par%dryslp;
         end if
         if(abs(dzbx(i,j))>dzmax .and. sum(dry(max(i-5,1):min(i+5,nx+1),j))>0.0d0) then    ! Robert search 5 cells left and right for dry point
            dzb=sign(1.0d0,dzbx(i,j))*(abs(dzbx(i,j))-dzmax)*(xz(i+1)-xz(i));;
            dzb=sign(1.0d0,dzb)*min(abs(dzb),0.005d0*par%dt*(xz(i+1)-xz(i))); !0.005d0 !Jaap, dzb_max = maximum allowed bottomchange per second per meter in direction of the avalanch 
            zb(i,j)=zb(i,j)+dzb*(xu(i+1)-xu(i))/(xu(i)-xu(i-1)); !Jaap make sure there is continuity of sediment in non uniform grids;
            sedero(i,j) = sedero(i,j)+dzb*(xu(i+1)-xu(i))/(xu(i)-xu(i-1));
            zs(i,j)=zs(i,j)+dzb*(xu(i+1)-xu(i))/(xu(i)-xu(i-1));
            zb(i+1,j)=zb(i+1,j)-dzb
            sedero(i+1,j)=sedero(i+1,j)-dzb
            zs(i+1,j)=zs(i+1,j)-dzb
         end if 
      end do
   end do
   dry = abs(1-wetu)
   do j=2,ny-1;
        do i=1,nx+1
           if(max(hh(i,j),hh(i,j+1))>par%hswitch+par%eps) then
              dzmax=par%wetslp;
           else
              dzmax=par%dryslp;
           end if
           if(abs(dzby(i,j))>dzmax .and. sum(dry(i,max(j-5,1):min(j+5,ny+1)))>0.0d0) then  ! Robert search 5 cells left and right for dry point
              dzb=sign(1.0d0,dzby(i,j))*(abs(dzby(i,j))-dzmax)*(yz(j+1)-yz(j));       
              dzb=sign(1.0d0,dzb)*min(abs(dzb),0.005d0*par%dt*(yz(j+1)-yz(j))); !0.005d0
              zb(i,j)=zb(i,j)+dzb*(yv(j+1)-yv(j))/(yv(j)-yv(j-1)); !Jaap make sure there is continuity of sediment in non uniform grids;;
              sedero(i,j)=sedero(i,j)+dzb*(yv(j+1)-yv(j))/(yv(j)-yv(j-1));
              zs(i,j)=zs(i,j)+dzb*(yv(j+1)-yv(j))/(yv(j)-yv(j-1));
              zb(i,j+1)=zb(i,j+1)-dzb
              sedero(i,j+1)=sedero(i,j+1)-dzb
              zs(i,j+1)=zs(i,j+1)-dzb
            end if
        end do
   end do
end do

endif

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine sb_vr(s,par)
use params
use spaceparams

IMPLICIT NONE

type(spacepars),target                  :: s
type(parameters)                        :: par

integer                                 :: i
integer                                 :: j,jg
real*8                                  :: dster,twothird
real*8                                  :: z0,Ass,cmax
real*8                                  :: Te,kvis,Sster,c1,c2,wster
real*8 , dimension(:,:),allocatable,save     :: vmag2,Cd,Asb,dhdx,dhdy,Ts
real*8 , dimension(:,:),allocatable,save     :: urms2,Ucr,term1,term2
real*8 , dimension(:,:),allocatable,save     :: uandv,b,fslope,hloc,kb,ceq

include 's.ind'
include 's.inp'

if (.not. allocated(vmag2)) then
   allocate (vmag2 (nx+1,ny+1))
   allocate (Cd    (nx+1,ny+1))
   allocate (Asb   (nx+1,ny+1))
   allocate (dhdx  (nx+1,ny+1))
   allocate (dhdy  (nx+1,ny+1))
   allocate (urms2 (nx+1,ny+1))
   allocate (Ucr   (nx+1,ny+1))
   allocate (term1 (nx+1,ny+1))
   allocate (term2 (nx+1,ny+1))
   allocate (uandv (nx+1,ny+1))
   allocate (b     (nx+1,ny+1))
   allocate (fslope(nx+1,ny+1))
   allocate (hloc  (nx+1,ny+1))
   allocate (kb    (nx+1,ny+1))
   allocate (Ts    (nx+1,ny+1))
   allocate (ceq    (nx+1,ny+1))
endif
! Soulsby van Rijn sediment transport formula
! Ad Reniers april 2006
!
! z is defined positive upward 
! x is defined positive toward the shore
!  Formal parameters:
!  ------------------
!
!   Var. I/O  Type Dimensions
!   -------------------------
!
! if (par%t == 150) pause
! cjaap: replaced par%hmin by par%eps
hloc   = max(hh,par%hmin) !Jaap par%hmin instead of par%eps
twothird=2.d0/3.d0
! use eulerian velocities
! cjaap: add turbulence near bottom
do j=1,ny+1 
    do i=1,nx+1
       kb(i,j) = (DR(i,j)/par%rho)**twothird/(exp( min( hloc(i,j)/max(H(i,j),0.1d0) ,100.)) -1.d0)
       vmag2(i,j) = ue(i,j)**2+ve(i,j)**2
    enddo
enddo

vmag2  = ue**2+ve**2
urms2  = urms**2+0.5*kb

do jg = 1,par%ngd

   ! cjaap: compute fall velocity with simple expression from Ahrens (2000)
   Te    = 20.d0
   kvis  = 4.d0/(20.d0+Te)*1d-5 ! Van rijn, 1993 
   Sster = s%D50(jg)/(4*kvis)*sqrt((par%rhos/par%rho-1)*par%g*s%D50(jg))
   c1    = 1.06d0*tanh(0.064d0*Sster*exp(-7.5d0/Sster**2))
   c2    = 0.22d0*tanh(2.34d0*Sster**-1.18d0*exp(-0.0064d0*Sster**2))
   wster = c1+c2*Sster
   par%w = wster*sqrt((par%rhos/par%rho-1.d0)*par%g*s%D50(jg))

   Ts       = 0.1d0*hh/par%w
   !Ts       = max(Ts,0.2d0)
   Tsg(:,:,jg) = max(Ts,0.2d0) 

   dster=25296*s%D50(jg)
   
   ! calculate treshold velocity Ucr
   if(s%D50(jg)<=0.0005) then
     Ucr=0.19*s%D50(jg)**0.1*log10(4*hloc/s%D90(jg))
   else if(s%D50(jg)<0.002) then
     Ucr=8.5*s%D50(jg)**0.6*log10(4*hloc/s%D90(jg))
   else
     write(*,*) ' s%D50(jg) > 2mm, out of validity range'
   end if
   ! drag coefficient
   z0 = par%z0
   Cd=(0.40/(log(max(hh,par%hmin)/z0)-1.0))**2 !Jaap
   !Cd = par%g/par%C**2;   ! consistent with flow modelling 

   ! transport parameters
   Asb=0.005*hloc*(s%D50(jg)/hloc/(1.65*par%g*s%D50(jg)))**1.2           ! bed load coefficent
   Ass=0.012*s%D50(jg)*dster**(-0.6)/(1.65*par%g*s%D50(jg))**1.2         ! suspended load coeffient
   term1=(vmag2+0.018/Cd*urms2)**0.5                                     ! nearbed-velocity
      
   term2 = 0*term1
   do j=1,ny+1
      do i=1,nx
         if(term1(i,j)>Ucr(i,j) .and. hh(i,j)>par%eps) then
            term2(i,j)=(term1(i,j)-Ucr(i,j))**2.4
!           term2(i,j)=(term1(i,j)-Ucr(i,j))**(1.7-0.7*tanh((ue(i,j)/sqrt(par%g*hh(i,j))-0.5)*10))
         end if
      end do
   end do
   ceq =(Asb+Ass)*term2  
   ceq = min(ceq,0.2)       ! equilibrium concentration
   ceq = ceq/hloc
   ceqg(:,:,jg) = ceq*sedcal(jg)
enddo  ! end og grain size classes

end subroutine