!!  Copyright (C)  Stichting Deltares, 2012-2014.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  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 registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine hdisp  ( pmsa   , fl     , ipoint , increm , noseg  ,
     &                    noflux , iexpnt , iknmrk , noq1   , noq2   ,
     &                    noq3   , noq4   )
!>\file
!>       (1D) Horizontal dispersion as velocity dependent reprofunction

C**********************************************************************
C     +----------------------------------------+
C     |    D E L F T   H Y D R A U L I C S     |
C     |    WAter Resources and Environment     |
C     +----------------------------------------+
C
C***********************************************************************
C
C     Project : OSIRIS 2002
C     Author  : Annette Kuin
C     Date    : 021216             Version : 0.01
C
C     History :
C
C     Date    Author          Description
C     ------  --------------  -----------------------------------
C     021216  Annette Kuin    Variable dispersion in 1D models
C
C***********************************************************************
C
C     Description of the module :
C
C Name    T   L I/O   Description                                    Units
C ----    --- -  -    -------------------                            -----

C     Logical Units : -

C     Modules called : -

C     Name     Type   Library
C     ------   -----  ------------

      IMPLICIT NONE
      REAL     PMSA  ( * ) , FL    (*)
      INTEGER  IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4

      INTEGER  IP1, IP2, IP3, IP4, IP5, IP6, IP7
      INTEGER  IN1, IN2, IN3, IN4, IN5, IN6, IN7
      REAL     VELOC, CHEZY, WIDTH, TOTDEP, alfaK,
     J         VELOCV, CHEZYV, WIDTHV, TOTDPV, alfaKV,
     J         VELOCN, CHEZYN, WIDTHN, TOTDPN, alfaKN,
     J         term1, term2, g, DVAR, MAXDSP, MAXDSN, MAXDSV
      INTEGER  IVAN  , INAAR , IQ

      IP1   = IPOINT( 1)
      IP2   = IPOINT( 2)
      IP3   = IPOINT( 3)
      IP4   = IPOINT( 4)
      IP5   = IPOINT( 5)
      IP6   = IPOINT( 6)
      IP7   = IPOINT( 7)
C
      IN1   = INCREM( 1)
      IN2   = INCREM( 2)
      IN3   = INCREM( 3)
      IN4   = INCREM( 4)
      IN5   = INCREM( 5)
      IN6   = INCREM( 6)
      IN7   = INCREM( 7)

c.....Exchangeloop over de verticale en breedterichtingen om
c.....ze op 0 te zetten en over de verticale richting om deze
c.....te initialiseren
      DO IQ=1, NOQ1+NOQ2+NOQ3+NOQ4
         PMSA(IP7 ) = 0.0
         IP7  = IP7  + IN7
      ENDDO

c.....Exchangeloop over horizontale lengterichting
      IP7 = IPOINT( 7)

      DO IQ = 1 , NOQ1

         IVAN  = IEXPNT(1,IQ)
         INAAR = IEXPNT(2,IQ)

         IF (IVAN.GT.0.OR.INAAR.GT.0) THEN

            IF (IVAN .LE. 0 ) IVAN = INAAR
            IF (INAAR .LE. 0 ) INAAR = IVAN

            VELOCV = PMSA(IP1 + (IVAN - 1) * IN1)
            WIDTHV = PMSA(IP2 + (IVAN - 1) * IN2)
            CHEZYV = PMSA(IP3 + (IVAN - 1) * IN3)
            TOTDPV = PMSA(IP4 + (IVAN - 1) * IN4)
            alfaKV = PMSA(IP5 + (IVAN - 1) * IN5)
            MAXDSV = PMSA(IP6 + (IVAN - 1) * IN6)

            VELOCN = PMSA(IP1 + (INAAR - 1) * IN1)
            WIDTHN = PMSA(IP2 + (INAAR - 1) * IN2)
            CHEZYN = PMSA(IP3 + (INAAR - 1) * IN3)
            TOTDPN = PMSA(IP4 + (INAAR - 1) * IN4)
            alfaKN = PMSA(IP5 + (INAAR - 1) * IN5)
            MAXDSN = PMSA(IP6 + (INAAR - 1) * IN6)

            VELOC  = (VELOCV + VELOCN) / 2.
            WIDTH  = (WIDTHV + WIDTHN) / 2.
            CHEZY  = (CHEZYV + CHEZYN) / 2.
            TOTDEP = (TOTDPV + TOTDPN) / 2.
            alfaK  = (alfaKV + alfaKN) / 2.
            MAXDSP = (MAXDSV + MAXDSN) / 2.

            g     = 9.81
            term1 = VELOC * WIDTH ** 2 * CHEZY
            term2 = TOTDEP * g ** 0.5
            DVAR = alfaK * term1 / term2

            !
            ! Limit the horizontal dispersion, if the value of
            ! MAXDSP on at least one side is positive
            !
            IF ( MAXDSV > 0.0 .AND. MAXDSN > 0.0 ) THEN
                PMSA(IP7) = MIN( DVAR, MAXDSP )
            ELSE
                MAXDSP = MAX( MAXDSV, MAXDSN )
                IF ( MAXDSP > 0.0 ) THEN
                    PMSA(IP7) = MIN( DVAR, MAXDSP )
                ELSE
                    PMSA(IP7) = DVAR
                ENDIF
            ENDIF

         ENDIF

         IP7  = IP7  + IN7

      ENDDO

      RETURN
      END