subroutine difuvl(icreep    ,timest    ,lundia    ,nst       ,icx       , &
                & icy       ,j         ,nmmaxj    ,nmmax     ,kmax      , &
                & lstsci    ,lstsc     ,lsal      ,ltem      ,lsecfl    , &
                & lsec      ,lsed      ,lsts      ,norow     ,irocol    , &
                & kcs       ,kcu       ,kfs       ,kfu       ,kfv       , &
                & kadu      ,kadv      ,s0        ,s1        ,hu        , &
                & hv        ,dps       ,qxk       ,qyk       ,qzk       , &
                & guu       ,gvv       ,guv       ,gvu       ,gsqs      , &
                & rbnd      ,sigdif    ,sigmol    ,r0        ,r1        , &
                & sour      ,sink      ,ws        ,thick     ,sig       , &
                & dicuv     ,vicww     ,dsdksi    ,dsdeta    ,dtdksi    , &
                & dtdeta    ,aak       ,bbk       ,cck       ,bdddx     , &
                & bddx      ,bdx       ,bux       ,buux      ,buuux     , &
                & uvdwk     ,vvdwk     ,areau     ,areav     ,volum0    , &
                & volum1    ,aakl      ,bbkl      ,cckl      ,ddkl      , &
                & bruvai    ,stage     ,eqmbcsand ,eqmbcmud  ,sedtyp    , &
                & seddif    ,gdp       )
!----- GPL ---------------------------------------------------------------------
!                                                                               
!  Copyright (C)  Stichting Deltares, 2011-2016.                                
!                                                                               
!  This program is free software: you can redistribute it and/or modify         
!  it under the terms of the GNU General Public License as published by         
!  the Free Software Foundation version 3.                                      
!                                                                               
!  This program 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 General Public License for more details.                                 
!                                                                               
!  You should have received a copy of the GNU General Public License            
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.        
!                                                                               
!  contact: delft3d.support@deltares.nl                                         
!  Stichting Deltares                                                           
!  P.O. Box 177                                                                 
!  2600 MH Delft, The Netherlands                                               
!                                                                               
!  All indications and logos of, and references to, "Delft3D" and "Deltares"    
!  are registered trademarks of Stichting Deltares, and remain the property of  
!  Stichting Deltares. All rights reserved.                                     
!                                                                               
!-------------------------------------------------------------------------------
!  $Id: difuvl.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/compute/difuvl.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Computes transport in the u, v and w-direction.
!              Explicit in the x- and y-direction, implicit in
!              sigma-direction.
!              Sinks are treated implicitly and sources explicit-
!              y. A special approach is used for the hori-
!              ontal diffusion to avoid artificial creeping.
! Method used: - Horizontal Advection in x- and y-direction :
!                explicit, Van Leer 2 scheme
!              - Horizontal Diffusion :
!                3D: explicit, along Z-planes
!                2D: explicit in U-direction
!                    explicit in V-direction
!              - Option: horizontal diffusion strictly horizontal
!                using special filter
!              - Vertical Advection :
!                implicit, central scheme
!              - Vertical Diffusion : implicit
!              - Sources are integrated explicitly.
!              - Sinks are integrated implicitly.
!              Reference: (B. Van Leer, Towards the Ultimate
!              Conservative Difference Scheme, II Monotonicity
!              and Conservation, combined in a second-order scheme,
!              J. Computational Phys. 14, 1974, 361-370)
!     Comment: For the Thatcher Harleman boundaries the boundary
!              points for outflow are reflected from the inner
!              points; for inflow the boundary conditions are
!              used (see also thahbc.for).
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use mathconsts
    use dfparall 
    use flow2d3d_timers
    use sediment_basics_module
   !
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    include 'flow_steps_f.inc'
    real(fp)                            , pointer :: ck
    real(fp)                            , pointer :: dicoww
    real(fp)                            , pointer :: eps
    real(fp), dimension(:,:,:)          , pointer :: fluxu
    real(fp), dimension(:,:,:)          , pointer :: fluxuc
    real(fp), dimension(:,:,:)          , pointer :: fluxv
    real(fp), dimension(:,:,:)          , pointer :: fluxvc
    real(fp)                            , pointer :: vicmol
    real(fp)                            , pointer :: xlo
    integer                             , pointer :: iro
    type (flwoutputtype)                , pointer :: flwoutput
!
! Global variables
!
    integer                                                     , intent(in)  :: icreep     !  Description and declaration in tricom.igs
    integer                                                     , intent(in)  :: icx        !  Increment in the X-dir., if ICX= NMAX then computation proceeds in the X-
                                                                                            !  dir. If icx=1 then computation proceeds in the Y-dir.
    integer                                                     , intent(in)  :: icy        !  Increment in the Y-dir. (see ICX)
    integer                                                     , intent(in)  :: j          !  Begin pointer for arrays which have been transformed into 1D arrays.
                                                                                            !  Due to the shift in the 2nd (M-)index, J = -2*NMAX + 1
    integer                                                     , intent(in)  :: kmax       !  Description and declaration in esm_alloc_int.f90
    integer                                                     , intent(in)  :: lsal       !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: lsec       !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: lsecfl     !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: lsed       !  Description and declaration in esm_alloc_int.f90
    integer                                                     , intent(in)  :: lsts       !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: lstsc      !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: lstsci     !  Description and declaration in esm_alloc_int.f90
    integer                                                     , intent(in)  :: ltem       !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: lundia     !  Description and declaration in inout.igs
    integer                                                     , intent(in)  :: nmmax      !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: nmmaxj     !  Description and declaration in dimens.igs
    integer                                                     , intent(in)  :: norow      !  Description and declaration in esm_alloc_int.f90
    integer                                                     , intent(in)  :: nst
    integer, dimension(5, norow)                                , intent(in)  :: irocol     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)                   , intent(in)  :: kcs        !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)                   , intent(in)  :: kcu        !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)                   , intent(in)  :: kfs        !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)                   , intent(in)  :: kfu        !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)                   , intent(in)  :: kfv        !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)             , intent(in)  :: kadu       !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)             , intent(in)  :: kadv       !  Description and declaration in esm_alloc_int.f90
    integer , dimension(lsed)                                   , intent(in)  :: sedtyp     !  sediment type: 0=total/1=noncoh/2=coh
    logical                                                     , intent(in)  :: eqmbcsand  !  Description and declaration in morpar.igs
    logical                                                     , intent(in)  :: eqmbcmud   !  Description and declaration in morpar.igs
    real(fp)                                                    , intent(in)  :: timest     !  Half Integration time step [sec.]
    real(prec), dimension(gdp%d%nmlb:gdp%d%nmub)                , intent(in)  :: dps        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: gsqs       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: guu        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: guv        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: gvu        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: gvv        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: hu         !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: hv         !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: s0         !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)                  , intent(in)  :: s1         !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)          , intent(in)  :: bruvai     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)          , intent(in)  :: vicww      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)          , intent(in)  :: qzk        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax, lsed)    , intent(in)  :: seddif     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax, lsed)    , intent(in)  :: ws         !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: aak        !  Internal work array (in CUCNP & UZD)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: bbk        !  Internal work array (in CUCNP & UZD)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: bdddx      !  Internal work array, implicit coupling of concentration in (N,M,K)
                                                                                            !  with layer concentration in (N,M-3,K)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: bddx       !  Internal work array, implicit coup-ing of concentration in (N,M,K)
                                                                                            !  with layer concentration in (N,M-2,K)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: bdx        !  Internal work array, implicit coupling of concentration in (N,M,K)
                                                                                            !  with layer concentration in (N,M-1,K)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: buuux      !  Internal work array, implicit coupling of concentration in (N,M,K)
                                                                                            !  with layer concentration in (N,M+3,K)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: buux       !  Internal work array, implicit coupling of concentration in (N,M,K)
                                                                                            !  with layer concentration in (N,M+2,K)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: bux        !  Internal work array, implicit coupling of concentration in (N,M,K)
                                                                                            !  with layer concentration in (N,M+1,K)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: cck        !  Internal work array (in CUCNP & UZD)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax+2)          , intent(in)  :: dicuv      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: dsdeta     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: dsdksi     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: dtdeta     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: dtdksi     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: areau
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: areav
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: qxk        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: qyk        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: volum0
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)            , intent(in)  :: volum1
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: uvdwk      !  Internal work array for Jac.iteration
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)                          :: vvdwk      !  Internal work array for Jac.iteration
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)                  :: aakl       !  Internal work array, lower diagonal tridiagonal matrix, implicit coupling
                                                                                            !  of concentration in (N,M,K) with concentration in (N,M,K-1)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)                  :: bbkl       !  Internal work array, main diagonal tridiagonal matrix, implicit coupling
                                                                                            !  of concentration in (N,M,K)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)                  :: cckl       !  Internal work array, upper diagonal tridiagonal matrix, implicit coupling
                                                                                            !  of concentration in (N,M,K) with concentration in (N,M,K+1)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)                  :: ddkl       !  Internal work array, diagonal space at (N,M,K,L)
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)                  :: r0         !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)                  :: r1         !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)    , intent(in)  :: sink       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci)    , intent(in)  :: sour       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)                                   , intent(in)  :: sig        !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)                                   , intent(in)  :: thick      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax, max(lstsc, 1), 2, norow)          , intent(in)  :: rbnd       !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(lstsci)                                 , intent(in)  :: sigdif     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(lstsci)                                 , intent(in)  :: sigmol     !  Description and declaration in esm_alloc_real.f90
    character(8)                                                , intent(in)  :: stage      !  First or second half time step
!
! Local variables
!
    integer                                                         :: ddb
    integer                                                         :: ic
    integer                                                         :: icxy
    integer                                                         :: istat
    integer                                                         :: k
    integer                                                         :: kfw
    integer                                                         :: l
    integer                                                         :: ll
    integer                                                         :: ls
    integer                                                         :: lst
    integer                                                         :: maskval
    integer                                                         :: mf
    integer                                                         :: ml
    integer                                                         :: n
    integer                                                         :: ndm
    integer                                                         :: nhystp
    integer                                                         :: nm
    integer                                                         :: nmd
    integer                                                         :: nmf
    integer                                                         :: nmfu
    integer                                                         :: nml
    integer                                                         :: nmlu
    integer                                                         :: nmu
    integer                                                         :: nmuu
    integer                                                         :: num
    integer                                                         :: nuum
    real(fp)                                                        :: adza
    real(fp)                                                        :: adzc
    real(fp)                                                        :: bi
    real(fp)                                                        :: cfl
    real(fp)                                                        :: cl
    real(fp)                                                        :: cr
    real(fp)                                                        :: ddzc
    real(fp)                                                        :: difiwe
    real(fp)                                                        :: difl
    real(fp)                                                        :: difr
    real(fp)                                                        :: diz1
    real(fp)                                                        :: flux
    real(fp)                                                        :: h0
    real(fp)                                                        :: h0i
    real(fp)                                                        :: h0new
    real(fp)                                                        :: h0old
    real(fp)                                                        :: qzw
    real(fp)                                                        :: r00
    real(fp), external                                              :: reddic
    real(fp)                                                        :: rr1
    real(fp)                                                        :: rr2
    real(fp)                                                        :: sqrtbv
    real(fp)                                                        :: tsg
    real(fp)                                                        :: uuu
    real(fp)                                                        :: vvv
    integer                                                         :: nm_pos ! indicating the array to be exchanged has nm index at the 2nd place, e.g., dbodsd(lsedtot,nm)
!
!! executable statements -------------------------------------------------------
!
    eps         => gdp%gdconst%eps
    fluxu       => gdp%gdflwpar%fluxu
    fluxuc      => gdp%gdflwpar%fluxuc
    fluxv       => gdp%gdflwpar%fluxv
    fluxvc      => gdp%gdflwpar%fluxvc
    flwoutput   => gdp%gdflwpar%flwoutput
    vicmol      => gdp%gdphysco%vicmol
    dicoww      => gdp%gdphysco%dicoww
    iro         => gdp%gdphysco%iro
    xlo         => gdp%gdturcoe%xlo
    ck          => gdp%gdturcoe%ck
    !
    ddb = gdp%d%ddbound
    icxy = max(icx, icy) 
    !
    istat = 0
    nm_pos= 1
    if (.not. associated(gdp%gdflwpar%fluxu)) then
        if (istat==0) allocate (gdp%gdflwpar%fluxu(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
        if (istat==0) allocate (gdp%gdflwpar%fluxv(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
        if (flwoutput%cumdifuflux) then
            if (istat==0) allocate (gdp%gdflwpar%fluxuc(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
            if (istat==0) allocate (gdp%gdflwpar%fluxvc(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
        endif
        !
        if (istat /= 0) then
            call prterr(lundia, 'U021', 'DIFUFVL: memory alloc error')
            call d3stop(1, gdp)
        endif
        !
        ! update local pointers
        !
        fluxu          => gdp%gdflwpar%fluxu
        fluxuc         => gdp%gdflwpar%fluxuc
        fluxv          => gdp%gdflwpar%fluxv
        fluxvc         => gdp%gdflwpar%fluxvc
        flwoutput      => gdp%gdflwpar%flwoutput
        !
        if (flwoutput%cumdifuflux) then
            fluxuc = 0.0
            fluxvc = 0.0
        endif
    endif
    !
    ! Initialization
    !
    fluxu = 0.0
    fluxv = 0.0
    !
    !  INITIALIZE FOR ALL (nm, k)
    !
    aak = 0.0
    cck = 0.0
    !
    do k = 1, kmax
       do nm = 1, nmmax
          if (kfs(nm)==1) then
             bbk(nm, k) = volum1(nm, k)/timest
          else
             bbk(nm, k) = 1.0
             if (lsec>0) r0(nm, k, lsecfl) = 0.0
          endif
       enddo
    enddo
    do l = 1, lstsci
       if (lsec==2 .and. l==lsecfl) then
          cycle
       endif
       do k = 1, kmax
          do nm = 1, nmmax
             if ( (kfs(nm)==1) .and. (kcs(nm)==1) ) then 
                ddkl(nm, k, l) = volum0(nm, k)*r0(nm, k, l)/timest
             else
                ddkl(nm, k, l) = r0(nm, k, l)
             endif
          enddo
       enddo
    enddo
    !
    nmd  = -icx
    nmu  = +icx
    nmuu = +icx + icx
    do nm = 1, nmmax
       !
       ! CONTRIBUTION OF ADVECTION IN X-DIRECTION
       !
       nmd  = nmd  + 1
       nmu  = nmu  + 1
       nmuu = nmuu + 1
       do k = 1, kmax
          if (kfu(nm)*kadu(nm, k) == 1) then
             uuu = qxk(nm, k)/(guu(nm)*hu(nm)*thick(k))
             cfl = uuu*timest/gvu(nm)
             if (qxk(nm, k) > 0.0) then
                do l = 1, lstsci
                   rr1 = abs(r0(nmd, k, l) - 2.0_fp*r0(nm, k, l) + r0(nmu, k, l))
                   rr2 = abs(r0(nmd, k, l) - r0(nmu, k, l))
                   if (kfu(nmd)*kadu(nmd, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp) then
                      r00 = r0(nm, k, l)
                   else
                      r00 = r0(nm, k, l)                                     &
                          & + (1.0_fp - cfl)*(r0(nm , k, l) - r0(nmd, k, l)) &
                          &                 *(r0(nmu, k, l) - r0(nm , k, l)) &
                          &                 /(r0(nmu, k, l) - r0(nmd, k, l))
                   endif
                   flux = qxk(nm, k)*r00
                   if (kcs(nm)  == 1) then                   
                        fluxu(nm,k,l) = fluxu(nm,k,l) + flux
                        ddkl (nm,k,l) = ddkl (nm,k,l) - flux
                   endif
                   if (kcs(nmu) == 1) ddkl(nmu, k, l) = ddkl(nmu, k, l) + flux
                enddo
             else
                do l = 1, lstsci
                   rr1 = abs(r0(nmuu, k, l) - 2.0_fp*r0(nmu, k, l) + r0(nm, k, l))
                   rr2 = abs(r0(nmuu, k, l) - r0(nm, k, l))
                   if (kfu(nmu)*kadu(nmu, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp) then
                      r00 = r0(nmu, k, l)
                   else
                      r00 = r0(nmu, k, l)                                     &
                          & + (1.0_fp + cfl)*(r0(nm , k, l) - r0(nmu , k, l)) &
                          &                 *(r0(nmu, k, l) - r0(nmuu, k, l)) &
                          &                 /(r0(nm , k, l) - r0(nmuu, k, l))
                   endif
                   flux = qxk(nm, k)*r00
                   if (kcs(nm)  == 1) then
                        fluxu(nm,k,l) = fluxu(nm,k,l) + flux
                        ddkl (nm,k,l) = ddkl (nm,k,l) - flux
                   endif
                   if (kcs(nmu) == 1) ddkl(nmu, k, l) = ddkl(nmu, k, l) + flux
                enddo
             endif
          endif
       enddo
    enddo
    !
    ! CONTRIBUTION OF ADVECTION IN Y-DIRECTION
    !
    ndm  = -icy
    num  = +icy
    nuum = +icy + icy
    do nm = 1, nmmax
       ndm  = ndm  + 1
       num  = num  + 1
       nuum = nuum + 1
       do k = 1, kmax
          if (kfv(nm)*kadv(nm, k) == 1) then
             vvv = qyk(nm, k)/(gvv(nm)*hv(nm)*thick(k))
             cfl = vvv*timest/guv(nm)
             if (qyk(nm, k) > 0.0) then
                do l = 1, lstsci
                   rr1 = abs(r0(ndm, k, l) - 2.0_fp*r0(nm, k, l) + r0(num, k, l))
                   rr2 = abs(r0(ndm, k, l) - r0(num, k, l))
                   if (kfv(ndm)*kadv(ndm, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp) then
                      r00 = r0(nm, k, l)
                   else
                      r00 = r0(nm, k, l)                                     &
                          & + (1.0_fp - cfl)*(r0(nm , k, l) - r0(ndm, k, l)) &
                          &                 *(r0(num, k, l) - r0(nm , k, l)) &
                          &                 /(r0(num, k, l) - r0(ndm, k, l))
                   endif
                   flux = qyk(nm, k) * r00
                   if (kcs(nm)  == 1) then
                        fluxv(nm,k,l) = fluxv(nm,k,l) + flux 
                        ddkl (nm,k,l) = ddkl (nm,k,l) - flux
                   endif
                   if (kcs(num) == 1) ddkl(num, k, l) = ddkl(num, k, l) + flux
                enddo
             else
                do l = 1, lstsci
                   rr1 = abs(r0(nuum, k, l) - 2.0_fp*r0(num, k, l) + r0(nm, k, l))
                   rr2 = abs(r0(nuum, k, l) - r0(nm, k, l))
                   if (kfv(num)*kadv(num, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp) then
                      r00 = r0(num, k, l)
                   else
                      r00 = r0(num, k, l)                                     &
                          & + (1.0_fp + cfl)*(r0(nm , k, l) - r0(num , k, l)) &
                          &                 *(r0(num, k, l) - r0(nuum, k, l)) &
                          &                 /(r0(nm , k, l) - r0(nuum, k, l))
                   endif
                   flux = qyk(nm, k) * r00
                   if (kcs(nm)  == 1) then 
                        fluxv(nm,k,l) = fluxv(nm,k,l) + flux
                        ddkl (nm,k,l) = ddkl (nm,k,l) - flux
                   endif
                   if (kcs(num) == 1) ddkl(num, k, l) = ddkl(num, k, l) + flux
                enddo
             endif
          endif
       enddo
    enddo
    !
    !
    ! Explicit algoritm (call DIFACR) leads to extra stablity criterium
    ! dt <= (dx**2)/(2*dicuv)
    !
    ! This diffusion part (loop 410) is constituent independent.
    ! The value of SIGDIF(L) = 0.7 (see TKECOF) for all LSTSCI
    !
    if (icreep==0 .or. kmax==1) then
       !
       ! HORIZONTAL DIFFUSION IN X-DIRECTION ALONG SIGMA PLANES
       !
       do l = 1, lstsci
          if (lsec==2 .and. l==lsecfl) then
             cycle
          endif
          do k = 1, kmax
             !
             ! CONTRIBUTION TO VOLUMES NM AND NMU
             !
             nmu = icx
             do nm = 1, nmmax
                nmu = nmu + 1
                if (kfu(nm)*kadu(nm, k) /= 0) then
                   cl             = r0   (nm ,k,l)
                   difl           = dicuv(nm ,k)
                   cr             = r0   (nmu,k,l)
                   difr           = dicuv(nmu,k)
                   flux           = 0.5* (cr-cl) * (difl+difr) / (0.7*gvu(nm))
                   maskval        = max(0, 2 - abs(kcs(nm))) 
                   fluxu(nm,k,l)  = fluxu(nm ,k,l) + areau(nm,k)*flux*maskval
                   ddkl (nm,k,l)  = ddkl (nm ,k,l) + areau(nm,k)*flux*maskval
                   maskval        = max(0, 2 - abs(kcs(nmu))) 
                   ddkl (nmu,k,l) = ddkl (nmu,k,l) - areau(nm,k)*flux*maskval
                endif
             enddo
          enddo
       enddo
       !
       ! HORIZONTAL DIFFUSION IN Y-DIRECTION ALONG SIGMA PLANES
       !
       do l = 1, lstsci
          if (lsec==2 .and. l==lsecfl) then
             cycle
          endif
          do k = 1, kmax
             !
             ! CONTRIBUTION TO VOLUMES NM AND NUM
             !
             num = icy
             do nm = 1, nmmax
                num = num + 1
                if (kfv(nm)*kadv(nm, k) /= 0) then
                   cl             = r0   (nm ,k,l)
                   difl           = dicuv(nm ,k)
                   cr             = r0   (num,k,l)
                   difr           = dicuv(num,k)
                   flux           = 0.5 * (cr-cl) * (difl+difr) / (0.7*guv(nm))
                   maskval        = max(0, 2 - abs(kcs(nm))) 
                   fluxv(nm ,k,l) = fluxv(nm ,k,l) + areav(nm,k)*flux*maskval
                   ddkl (nm ,k,l) = ddkl (nm ,k,l) + areav(nm,k)*flux*maskval
                   maskval        = max(0, 2 - abs(kcs(num))) 
                   ddkl (num,k,l) = ddkl (num,k,l) - areav(nm,k)*flux*maskval
                endif
             enddo
          enddo
       enddo
    else
       !
       ! Explicit algoritm (call DIFACR) leads to extra stablity criterium
       ! dt <= (dx**2)/(2*dicuv)
       !
       ! HORIZONTAL DIFFUSION ALONG Z-PLANES (only if KMAX > 1 and Anti Creep)
       !
       call difacr(icx       ,icy       ,j         ,nmmaxj    ,nmmax     , &
                 & kmax      ,lstsci    ,lsal      ,ltem      ,kcs       , &
                 & kfu       ,kfv       ,kadu      ,kadv      ,s0        , &
                 & dps       ,r0        ,ddkl      ,guu       ,gvv       , &
                 & guv       ,gvu       ,thick     ,sig       ,dicuv     , &
                 & sigdif    ,dsdksi    ,dtdksi    ,dsdeta    ,dtdeta    , &
                 & gdp       )
    endif
    !
    ! Swap fluxes when in stage 2, because then the U- and V- direction are interchanged
    !
    if (stage == 'stage2  ') then
       fluxu       => gdp%gdflwpar%fluxv
       fluxv       => gdp%gdflwpar%fluxu
    endif
    !
    ! Cumulative flux
    !
    if (flwoutput%cumdifuflux) then
       fluxuc = fluxuc + fluxu * timest
       fluxvc = fluxvc + fluxv * timest
    endif
    !
    if (kmax>1) then
       do k = 1, kmax - 1
          if (k==1 .or. k==kmax - 1) then
             kfw = 1
          else
             kfw = 0
          endif
          do nm = 1, nmmax
             !
             ! ADVECTION IN VERTICAL DIRECTION; W*DU/DZ
             !
             if (kfs(nm)==1) then
                qzw = qzk(nm, k)
                if (qzw>0.0) then
                   adza = 0.5*qzw*(1 - kfw)
                   adzc = 0.5*qzw*(1 + kfw)
                else
                   adza = 0.5*qzw*(1 + kfw)
                   adzc = 0.5*qzw*(1 - kfw)
                endif
                aak(nm, k + 1) = aak(nm, k + 1) + adza
                bbk(nm, k + 1) = bbk(nm, k + 1) + adzc
                bbk(nm, k) = bbk(nm, k) - adza
                cck(nm, k) = cck(nm, k) - adzc
             endif
          enddo
       enddo
    endif
    do l = 1, lstsci
       if (lsec==2 .and. l==lsecfl) then
          cycle
       endif
       do k = 1, kmax
          do nm = 1, nmmax
             if (kfs(nm)==1) then
                aakl(nm, k, l) = aak(nm, k)
                bbkl(nm, k, l) = bbk(nm, k)
                cckl(nm, k, l) = cck(nm, k)
             endif
          enddo
       enddo
    enddo
    !
    if (kmax>1) then
       do l = 1, lstsci
          !
          ! l = sediment: ls > 0
          ! else        : ls = 0
          !
          ls = 0
          if ((l>max(lsal, ltem)) .and. (l<=lsts)) then
             ls = l - max(lsal, ltem)
          endif
          do k = 1, kmax - 1
             tsg = 0.5*(thick(k) + thick(k + 1))
             do nm = 1, nmmax
                !
                ! DIFFUSION IN VERTICAL DIRECTION
                !
                if (kfs(nm) == 1) then
                   h0  = max(0.1_fp, s0(nm) + real(dps(nm),fp))
                   h0i = 1.0 / h0
                   !
                   ! Internal wave contribution
                   !
                   sqrtbv = max(0.0_fp, bruvai(nm, k))
                   sqrtbv = sqrt(sqrtbv)
                   difiwe = 0.2 * sqrtbv * xlo**2
                   if (ls > 0) then
                      !
                      ! sediment constituent:
                      ! No dicoww-restriction in reddic
                      !
                      diz1 = vicmol/sigmol(l) + difiwe + seddif(nm, k, ls)/sigdif(l)
                   else
                      !
                      ! all other constituents:
                      ! dicoww-restriction is moved from TURCLO to here (in reddic)
                      ! vicww is used instead of dicww
                      !
                      diz1 = vicmol/sigmol(l) + reddic(difiwe + vicww(nm,k)/sigdif(l), gdp)
                   endif
                   ddzc = gsqs(nm) * diz1 * h0i / tsg
                   aakl(nm, k+1, l) = aakl(nm, k+1, l) - ddzc
                   bbkl(nm, k+1, l) = bbkl(nm, k+1, l) + ddzc
                   bbkl(nm, k  , l) = bbkl(nm, k  , l) + ddzc
                   cckl(nm, k  , l) = cckl(nm, k  , l) - ddzc
                endif
             enddo
          enddo
       enddo
    endif
    !
    ! Include settling velocities and Dirichlet BC for sediments in
    ! matrices AAKL/BBKL/CCKL/DDKL
    !
    if (lsed > 0) then
       call timer_start(timer_difu_difws, gdp)
       call dif_ws(j         ,nmmaxj    ,nmmax     ,kmax      ,lsal      , &
                 & ltem      ,lstsci    ,lsed      ,kcs       ,kfs       , &
                 & gsqs      ,ws        ,aakl      ,bbkl      ,cckl      , &
                 & gdp       )
       call timer_stop(timer_difu_difws, gdp)
    endif
    !
    ! SET VALUES IN OPEN BOUNDARY POINTS (IN PART. FOR Y-DIRECTION)
    !     On open boundary no secondary flow (=> loop over LSTSC)
    !
    do nm = 1, nmmax
       if (kcs(nm)==2) then
          do l = 1, lstsc
             do k = 1, kmax
                ddkl(nm, k, l) = r0(nm, k, l)
                aakl(nm, k, l) = 0.0
                bbkl(nm, k, l) = 1.0
                cckl(nm, k, l) = 0.0
             enddo
          enddo
       endif
    enddo
    !
    ! BOUNDARY CONDITIONS
    !     On open boundary no secondarY flow (=> loop over LSTSC)
    !
    do ic = 1, norow
       n = irocol(1, ic)
       mf = irocol(2, ic) - 1
       ml = irocol(3, ic)
       nmf = (n + ddb)*icy + (mf + ddb)*icx - icxy
       nml = (n + ddb)*icy + (ml + ddb)*icx - icxy
       nmfu = nmf + icx
       nmlu = nml + icx
       !
       ! IMPLEMENTATION OF BOUNDARY CONDITIONS
       !
       if (kcu(nmf)==1) then
          do k = 1, kmax
             do l = 1, lstsc
                ddkl(nmf, k, l) = rbnd(k, l, 1, ic)
             enddo
          enddo
       endif
       if (kcu(nml)==1) then
          do k = 1, kmax
             do l = 1, lstsc
                ddkl(nmlu, k, l) = rbnd(k, l, 2, ic)
             enddo
          enddo
       endif
       !
       ! optional Neumann boundary condition for suspended sediment fractions
       !
       lst = max(lsal, ltem)
       do l = 1, lsed
          ll = lst + l
          if ((eqmbcsand .and. sedtyp(l) == SEDTYP_NONCOHESIVE_SUSPENDED) .or. &
            & (eqmbcmud  .and. sedtyp(l) == SEDTYP_COHESIVE)             ) then
             if (kcu(nmf) == 1) then
                do k = 1, kmax
                   ddkl(nmf, k, ll) = max(0.0_fp, r0(nmfu, k, ll))
                enddo
             endif
             if (kcu(nml) == 1) then
                do k = 1, kmax
                   ddkl(nmlu, k, ll) = max(0.0_fp, r0(nml, k, ll))
                enddo
             endif
          endif
       enddo
    enddo
    do l = 1, lstsci
       !
       ! SOURCES AND SINK TERMS
       !
       ! SINKS ARE TREATED IMPLICITLY
       !
       if (lsec==2 .and. l==lsecfl) then
          !
          ! secondary flow (equilibrium equals to new intensity)
          !
          do k = 1, kmax
             do nm = 1, nmmax
                if ( (kfs(nm)==1) .and. (kcs(nm)==1) ) then 
                   h0new = s1(nm) + real(dps(nm),fp)
                   if (abs(sink(nm, k, l))*h0new>eps) then
                      h0old = s0(nm) + real(dps(nm),fp)
                      r1(nm, k, l) = sour(nm, k, l)*h0old/(sink(nm, k, l)*h0new)
                   else
                      r1(nm, k, l) = 0.0
                   endif
                endif
             enddo
          enddo
       else
          do k = 1, kmax
             do nm = 1, nmmax
                if ( (kfs(nm)==1) .and. (kcs(nm)==1) ) then 
                   bbkl(nm, k, l) = bbkl(nm, k, l) + sink(nm, k, l)*volum1(nm, k)
                   ddkl(nm, k, l) = ddkl(nm, k, l) + sour(nm, k, l)*volum0(nm, k)
                endif
             enddo
          enddo
          !
          ! set concentraties in dry points and in open boundary points
          !
          do k = 1, kmax
             do nm = 1, nmmax
                if ((kfs(nm)==0 .and. kcs(nm)==1) .or. kcs(nm)==2) then
                   r1(nm, k, l) = ddkl(nm, k, l)
                endif
             enddo
          enddo
       endif
       !
       if (l==lsecfl) then
          !
          ! boundary conditions secondary flow (spiral motion intensity)
          !
          call secbou(j         ,nmmaxj    ,kmax      ,icx       ,icy       , &
                    & lstsci    ,lsecfl    ,kfu       ,irocol    ,norow     , &
                    & s0        ,s1        ,dps       ,r1        ,sour      , &
                    & sink      ,gdp       )
          if (lsec==2) then
             ! 
             ! exchange r1 with neighbours for parallel runs 
             ! 
             call dfexchg ( r1(:,:,l), 1, kmax, dfloat, nm_pos, gdp ) 
             cycle
          endif
       endif
       !
       ! DD code added:
       !
       ! left-hand side is now set by Delft3D-FLOW instead of the mapper
       !
       do nm = 1, nmmax
          if (kcs(nm) == 3 ) then
             do k = 1, kmax
                aakl(nm,k,l) = 0.0
                bbkl(nm,k,l) = 1.0
                cckl(nm,k,l) = 0.0
                ddkl(nm,k,l) = r0(nm,k,l)
             enddo
          endif
       enddo
       !
       !
       !        D3dFlow_Build_ADI_Conc: poke the coupling equations into system
       !
       nhystp = nxtstp(d3dflow_build_adi_conc, gdp)
       !
       ! DD code added end
       !
       !***SOLUTION PROCEDURE SYSTEM OF EQUATIONS
       !
       !
       ! DD code added:
       !
       !
       !       (re)solve system of equations
       !
  111  continue
       gdp%dd%difuiter = gdp%dd%difuiter + 1
       !
       ! set concentrations in coupling points
       !
       do k = 1, kmax
          do nm = 1, nmmax
             if (kcs(nm)==3) then
                r1(nm, k, l) = ddkl(nm, k, l)
             endif
          enddo
       enddo
       !
       ! DD code added end
       !
       do nm = 1, nmmax
          if ( (kfs(nm)==1) .and. (kcs(nm)==1) ) then 
             !
             !***SCALE ROWS OF MATRIX/RIGHT HAND SIDE VECTOR
             !
             do k = 1, kmax
                bi             = 1.0_fp / bbkl(nm, k, l)
                aakl(nm, k, l) = aakl(nm, k, l) * bi
                bbkl(nm, k, l) = 1.0_fp
                cckl(nm, k, l) = cckl(nm, k, l) * bi
                ddkl(nm, k, l) = ddkl(nm, k, l) * bi
             enddo
             r1(nm, 1, l) = ddkl(nm, 1, l)
             !
             !***SOLUTION PROCEDURE SYSTEM OF EQUATIONS
             !
             !
             ! Division by the pivot for k=1 is not needed anymore
             ! because of row scaling
             !
             do k = 2, kmax
                bi = 1./(bbkl(nm, k, l) - aakl(nm, k, l)*cckl(nm, k - 1, l))
                cckl(nm, k, l) = cckl(nm, k, l) * bi
                r1(nm, k, l) = (ddkl(nm, k, l) - aakl(nm, k, l)*r1(nm, k - 1, l)) * bi
             enddo
             do k = kmax - 1, 1, -1
                r1(nm, k, l) = r1(nm, k, l) - cckl(nm, k, l)*r1(nm, k + 1, l)
             enddo
          endif
       enddo
       !
       ! exchange r1 with neighbours for parallel runs 
       ! 
       call dfexchg ( r1(:,:,l), 1, kmax, dfloat, nm_pos, gdp ) 
       ! 
       ! DD code added:
       !
       !
       !       D3dFlow_Solve_ADI_Conc: Check for convergence
       !
       nhystp = nxtstp(d3dflow_solve_adi_conc, gdp)
       if (nhystp==d3dflow_solve_adi_conc) goto 111
       !
       ! DD code added end
       !
    enddo
    !
    ! Deallocate possible local array
    !
end subroutine difuvl
