text:
!!****h* Conquest/rng *
!!  NAME
!!   rng
!!  PURPOSE
!!   A place for all random number generator-related code. Mersenne Twister
!!   code adapted from original C code by Takuji Nishimura and Makoto
!!   Matsumoto, and translated to Fortran by Remi Piatek (see below)
!!  AUTHOR
!!   Zamaan Raza
!!  CREATION DATE
!!   2019/05/16
!!  MODIFICATION HISTORY
!!   2019/05/21 zamaan
!!    Changes to intialisation to make rng_normal easier to use. Now one call
!!    returns a single random number, (the other is buffered) but it needs to
!!    be initialised using init_normal
!!  SOURCE
!!
!-------------------------------------------------------------------------------
!   This is a Fortran translation of the 64-bit version of
!   the Mersenne Twister pseudorandom number generator
!
!   Translated from C-program for MT19937-64 (2004/9/29 version)
!   originally coded by Takuji Nishimura and Makoto Matsumoto
!   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html
!
!   Fortran translation by Remi Piatek
!   The University of Copenhagen
!   Department of Economics
!   email: {first}.{last}@econ.ku.dk
!
!-------------------------------------------------------------------------------
!   A C-program for MT19937-64 (2004/9/29 version).
!   Coded by Takuji Nishimura and Makoto Matsumoto.
!
!   This is a 64-bit version of Mersenne Twister pseudorandom number
!   generator.
!
!   Before using, initialize the state by using init_genrand64(seed)
!   or init_by_array64(init_key, key_length).
!
!   Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
!   All rights reserved.
!
!   Redistribution and use in source and binary forms, with or without
!   modification, are permitted provided that the following conditions
!   are met:
!
!     1. Redistributions of source code must retain the above copyright
!        notice, this list of conditions and the following disclaimer.
!
!     2. Redistributions in binary form must reproduce the above copyright
!        notice, this list of conditions and the following disclaimer in the
!        documentation and/or other materials provided with the distribution.
!
!     3. The names of its contributors may not be used to endorse or promote
!        products derived from this software without specific prior written
!        permission.
!
!   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
!   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
!   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
!   OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
!   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
!   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
!   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
!   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!   References:
!   T. Nishimura, ``Tables of 64-bit Mersenne Twisters''
!     ACM Transactions on Modeling and
!     Computer Simulation 10. (2000) 348--357.
!   M. Matsumoto and T. Nishimura,
!     ``Mersenne Twister: a 623-dimensionally equidistributed
!       uniform pseudorandom number generator''
!     ACM Transactions on Modeling and
!     Computer Simulation 8. (Jan. 1998) 3--30.
!
!   Any feedback is very welcome.
!   http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html
!   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)
!-------------------------------------------------------------------------------
module rng

  use datatypes

  implicit none

  integer(wide), parameter, private :: nn       = 312_wide
  integer(wide), parameter, private :: mm       = 156_wide
  integer(wide), parameter, private :: seed_def = 5489_wide
  integer(wide), parameter, private :: matrix_a = -5403634167711393303_wide
  integer(wide), parameter, private :: um       = -2147483648_wide
  ! most significant 33 bits
  integer(wide), parameter, private :: lm       = -2147483647_wide
  ! least significant 31 bits
  real(double), parameter, private      :: pi253_1 = 1._double/(2._double**53 -&
                                                     1._double)
  real(double), parameter, private      :: pi253   = 1._double/(2._double**53)
  real(double), parameter, private      :: pi252   = 1._double/(2._double**52)

  type type_rng

    integer(wide), allocatable  :: seed(:)

    integer(wide) :: mt(nn)     ! array for state vector
    integer       :: mti = nn+1 ! mti==nn+1 means mt(nn) is not initialized

    ! For normal distribution
    logical       :: buffered
    real(double)  :: buffer, sigma, mu

    contains

      procedure, private  :: init_genrand64
      procedure, private  :: init_by_array64
      procedure, private  :: genrand64_int64
      procedure, private  :: genrand64_real1
      procedure, private  :: genrand64_real2
      procedure, private  :: genrand64_real3
      procedure, private  :: get_random_seed
      procedure, public   :: init_rng
      procedure, public   :: init_normal
      procedure, public   :: rng_uniform
      procedure, public   :: rng_normal
      procedure, public   :: rng_integer
  end type type_rng
  !!***

  contains

    !!****f* rng/init_genrand64 *
    !!
    !!  NAME
    !!   init_genrand64
    !!  PURPOSE
    !!   Initialises mt(nn) with a seed
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    subroutine init_genrand64(rn, seed)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn
      integer(wide), intent(in)       :: seed

      ! Local variables
      integer                 :: i

      rn%mt(1) = seed
      do i = 1, nn-1
        rn%mt(i+1) = 6364136223846793005_wide * ieor(rn%mt(i), ishft(rn%mt(i), -62)) + i
      end do

      rn%mti = nn

    end subroutine init_genrand64
    !!***

    !!****f* rng/init_by_array64 *
    !!
    !!  NAME
    !!   init_by_array64
    !!  PURPOSE
    !!   Initialises by array with an array length
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    subroutine init_by_array64(rn, init_key)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn
      integer(wide), intent(in)       :: init_key(:)

      ! Local variables
      integer(wide), parameter  :: c1 = 3935559000370003845_wide
      integer(wide), parameter  :: c2 = 2862933555777941757_wide
      integer(wide)             :: i, j, k, kk, key_length

      call rn%init_genrand64(19650218_wide)
      key_length = size(init_key)
      i = 1_wide; j = 0_wide
      k = max(nn, key_length)

      do kk = 1, k
        rn%mt(i+1) = ieor(rn%mt(i+1), c1 * ieor(rn%mt(i), ishft(rn%mt(i), -62))) &
                    + init_key(j+1) + j
        i = i+1; j = j+1
        if(i >= nn) then
          rn%mt(1) = rn%mt(nn)
          i = 1
        end if
        if(j >= key_length) j = 0
      end do

      do kk = 1, nn-1
        rn%mt(i+1) = ieor(rn%mt(i+1), c2 * ieor(rn%mt(i), ishft(rn%mt(i), -62))) - i
        i = i+1
        if(i >= nn) then
          rn%mt(1) = rn%mt(nn)
          i = 1
        end if
      end do

      rn%mt(1) = ishft(1_wide, 63)  ! MSB is 1; assuring non-zero initial array

    end subroutine init_by_array64
    !!***

    !!****f* rng/genrand64_int64 *
    !!
    !!  NAME
    !!   genrand64_int64
    !!  PURPOSE
    !!   Generates random integer on the interval [-2^63,2^63-1]
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    integer(double) function genrand64_int64(rn)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn

      ! Local variables
      integer(wide) :: mag01(0:1) = (/0_wide, matrix_a/)
      integer(wide) :: x
      integer     :: i

      if(rn%mti >= nn) then ! generate nn words at one time

        ! if init_genrand64() has not been called, a default initial seed is used
        if(rn%mti == nn+1) call rn%init_genrand64(seed_def)

        do i = 1, nn-mm
          x = ior(iand(rn%mt(i),um), iand(rn%mt(i+1), lm))
          rn%mt(i) = ieor(ieor(rn%mt(i+mm), ishft(x, -1)), mag01(iand(x, 1_wide)))
        end do

        do i = nn-mm+1, nn-1
          x = ior(iand(rn%mt(i), um), iand(rn%mt(i+1), lm))
          rn%mt(i) = ieor(ieor(rn%mt(i+mm-nn), ishft(x, -1)), mag01(iand(x, 1_wide)))
        end do

        x = ior(iand(rn%mt(nn), um), iand(rn%mt(1), lm))
        rn%mt(nn) = ieor(ieor(rn%mt(mm), ishft(x, -1)), mag01(iand(x, 1_wide)))

        rn%mti = 0

      end if

      rn%mti = rn%mti + 1
      x = rn%mt(rn%mti)

      x = ieor(x, iand(ishft(x,-29), 6148914691236517205_wide))
      x = ieor(x, iand(ishft(x, 17), 8202884508482404352_wide))
      x = ieor(x, iand(ishft(x, 37),   -2270628950310912_wide))
      x = ieor(x, ishft(x, -43))

      genrand64_int64 = x

    end function genrand64_int64
    !!***

    !!****f* rng/genrand64_real1 *
    !!
    !!  NAME
    !!   genrand64_real1
    !!  PURPOSE
    !!   Generates a random real number on the interval [0,1]
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    real(double) function genrand64_real1(rn)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn

        genrand64_real1 = real(ishft(rn%genrand64_int64(), -11), kind=double) * pi253_1

    end function genrand64_real1

    !!****f* rng/genrand64_real2 *
    !!
    !!  NAME
    !!   genrand64_real2
    !!  PURPOSE
    !!   Generates a random real number on the interval [0,1)
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    real(double) function genrand64_real2(rn)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn

        genrand64_real2 = real(ishft(rn%genrand64_int64(), -11), kind=double) * pi253

    end function genrand64_real2
    !!***

    !!****f* rng/genrand64_real3 *
    !!
    !!  NAME
    !!   genrand64_real3
    !!  PURPOSE
    !!   Generates a random real number on the interval (0,1)
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    real(double) function genrand64_real3(rn)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn

      genrand64_real3 = real(ishft(rn%genrand64_int64(), -12), kind=double)
      genrand64_real3 = (genrand64_real3 + 0.5_double) * pi252

    end function genrand64_real3
    !!***

    !!****f* rng/get_random_seed *
    !!
    !!  NAME
    !!   get_random_seed
    !!  PURPOSE
    !!   Generates a random real number on the interval (0,1)
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!   2019/05/22 zamaan
    !!    Removed the compiler-dependent getpid() call
    !!  SOURCE
    !!
    subroutine get_random_seed(rn)

      use input_module, only: io_assign, io_close

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn

      ! Local variables
      integer               :: i, un, istat, dt(8), pid, t(2), s, n
      integer(8)            :: count, tms

      call random_seed(size=n)
      allocate(rn%seed(n))
      ! First try if the OS provides a random number generator
      call io_assign(un)
      open(unit=un, file="/dev/urandom", access="stream", &
           form="unformatted", action="read", status="old", iostat=istat)
      if (istat == 0) then
         read(un) rn%seed
         close(un)
      else
         ! Fallback to XOR:ing the current time and pid. The PID is
         ! useful in case one launches multiple instances of the same
         ! program in parallel.
         call system_clock(count)
         if (count /= 0) then
            t = transfer(count, t)
         else
            call date_and_time(values=dt)
            tms = (dt(1) - 1970) * 365_8 * 24 * 60 * 60 * 1000 &
                 + dt(2) * 31_8 * 24 * 60 * 60 * 1000 &
                 + dt(3) * 24 * 60 * 60 * 60 * 1000 &
                 + dt(5) * 60 * 60 * 1000 &
                 + dt(6) * 60 * 1000 + dt(7) * 1000 &
                 + dt(8)
            t = transfer(tms, t)
         end if
         s = ieor(t(1), t(2))
         ! getpid is compiler-dependent; intrinsic in gfortran but but requires
         ! ifport for ifort - zamaan
!         pid = getpid() + 1099279 ! Add a prime
!         s = ieor(s, pid)
         if (n >= 3) then
            rn%seed(1) = t(1) + 36269
            rn%seed(2) = t(2) + 72551
!            rn%seed(3) = pid
!            If (n > 3) then
!               rn%seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /)
            if (n > 2) then
               rn%seed(3:) = s + 37 * (/ (i, i = 0, n - 3) /)
            end if
         else
            rn%seed = s + 37 * (/ (i, i = 0, n - 1 ) /)
         end if
      end if
      call io_close(un)

    end subroutine get_random_seed
    !!***

    !!****f* rng/init_rng *
    !!
    !!  NAME
    !!   init_rng
    !!  PURPOSE
    !!   Conquest wrapper for seed generator and MT rng. Initialise using
    !!   int_seed if provided, otherwise use get_random_seed
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!   2019/05/21 zamaan
    !!    Now takes rng_seed from global_module, generates a random seed array
    !     if rng_seed < 0 (default = -1)
    !!  SOURCE
    !!
    subroutine init_rng(rn)

      use global_module, only: rng_seed

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)    :: rn

      if (rng_seed > 0)  then
        call init_genrand64(rn, int(rng_seed,wide))
      else
        call rn%get_random_seed
        call rn%init_by_array64(rn%seed)
      end if
    end subroutine init_rng
    !!***

    !!****f* rng/init_normal *
    !!
    !!  NAME
    !!   init_normal
    !!  PURPOSE
    !!   Initialise generator for normal distribution
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/21
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    subroutine init_normal(rn, sigma, mu)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn
      real(double), intent(in)        :: sigma, mu

      rn%buffered = .false.
      rn%sigma = sigma
      rn%mu = mu

    end subroutine init_normal
    !!***

    !!****f* rng/rng_uniform *
    !!
    !!  NAME
    !!   rng_uniform
    !!  PURPOSE
    !!   Conquest wrapper for real random number generator.
    !!  USAGE
    !!   interval 1 (default): [0,1], interval 2: [0,1), interval 3: (0,1)
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/16
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    real(double) function rng_uniform(rn, interval)

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)    :: rn
      integer, optional, intent(in)     :: interval

      ! Local variables
      real(double)                      :: num

      if (present(interval)) then
        select case(interval)
        case(1)
          num = rn%genrand64_real1()
        case(2)
          num = rn%genrand64_real2()
        case(3)
          num = rn%genrand64_real3()
        end select
      else
        num = rn%genrand64_real1()
      end if

      rng_uniform = real(num, double)

    end function rng_uniform
    !!***

    !!****f* move_atoms/rng_normal *
    !!  NAME
    !!   rng_normal
    !!  PURPOSE
    !!   Box-Muller transform to generate pairs of numbers drawn from a
    !!   Gaussian distribution given a uniform distribution, via
    !!   Marsaglia's polar method (no sine or cosine calls)
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/17
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    real(double) function rng_normal(rn)

      use numbers

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn

      ! Local variables
      real(double)              :: x1, x2, y1, y2, w

      if (rn%buffered) then
        rng_normal = rn%buffer
        rn%buffered = .false.
      else
        do
          x1 = rn%rng_uniform()
          x2 = rn%rng_uniform()
          y1 = two*x1 - one
          y2 = two*x2 - one
          w = y1*y1 + y2*y2
          if (w < one) exit
        end do
        w = sqrt(-two*log(w)/w)
        rng_normal = y1*w*rn%sigma + rn%mu
        rn%buffer  = y2*w*rn%sigma + rn%mu
        rn%buffered = .true.
      end if

    end function rng_normal
    !!***

    !!****f* move_atoms/rng_integer *
    !!  NAME
    !!   rng_integer
    !!  PURPOSE
    !!   Generate random numbers in the interval [a,b)
    !!  AUTHOR
    !!   Zamaan Raza
    !!  CREATION DATE
    !!   2019/05/17
    !!  MODIFICATION HISTORY
    !!
    !!  SOURCE
    !!
    integer function rng_integer(rn, a, b)

      use datatypes

      implicit none

      ! Passed variables
      class(type_rng), intent(inout)  :: rn
      integer, intent(in)             :: a, b

      ! Local variables
      integer :: rint

      rint = modulo(rn%genrand64_int64(), (b-a))
      rng_integer = int(rint) + a

    end function  rng_integer
    !!***

end module rng
!!***

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
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

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd