bedload_hunz_meyer.f

Go to the documentation of this file.

!                   *****************************
                    SUBROUTINE bedload_hunz_meyer
!                   *****************************
!
     &  (tob, mu, acladm, unladm, npoin, dens, xmve, grav, dm, ac,
     &    tetap, ahunzi, acp, hiding, qsc)
!
!***********************************************************************
! SISYPHE   V6P1                                   21/07/2011
!***********************************************************************
!
!brief    HUNZIKER BEDLOAD FORMULATION (1995)
!+
!+           (ADAPTED FROM MEYER-PETER FORMULATION).
!
!note     **-1.5 AND **1.5 SHOULD BE OPTIMISED (JMH)
!
!history  BUI MINH DUC
!+        **/01/2002
!+        V5P2
!+
!
!history  C. VILLARET
!+        **/10/2003
!+        V5P4
!+
!
!history
!+
!+        V6P0
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  C.VILLARET (EDF-LNHE), P.TASSI (EDF-LNHE)
!+        19/07/2011
!+        V6P1
!+  Name of variables
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AC             |<->| CRITICAL SHIELDS PARAMETER
!| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
!| ACP            |<->| MODIFIED SHIELDS PARAMETER
!| AHUNZI         |<->| COEFFICIENT OF HUNZIKER FORMULA
!| DENS           |-->| RELATIVE DENSITY
!| DM             |-->| SEDIMENT GRAIN DIAMETER
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| HIDING         |-->| HIDING FACTOR CORRECTION
!| MU             |<->| CORRECTION FACTOR FOR BED ROUGHNESS
!| NPOIN          |-->| NUMBER OF POINTS
!| QSC            |<->| BED LOAD TRANSPORT
!| TETAP          |<->| DIMENSIONLESS BED SHEAR STRESS
!| TOB            |<->| BED SHEAR STRESS (TOTAL FRICTION)
!| UNLADM         |-->| MEAN DIAMETER OF ACTIVE STRATUM LAYER
!| XMVE           |-->| FLUID DENSITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,
     &          ex_bedload_hunz_meyer => bedload_hunz_meyer
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
      ! 2/ GLOBAL VARIABLES
      ! -------------------
      TYPE(bief_obj),   INTENT(IN)    :: TOB, MU, ACLADM, UNLADM
      INTEGER,          INTENT(IN)    :: NPOIN
      DOUBLE PRECISION, INTENT(IN)    :: DENS, XMVE, GRAV, DM, AC
      TYPE(bief_obj),   INTENT(INOUT) :: TETAP, AHUNZI ! WORK ARRAY T1, T2
      TYPE(bief_obj),   INTENT(INOUT) :: ACP           ! WORK ARRAY T3
      TYPE(bief_obj),   INTENT(INOUT)   :: HIDING, QSC
!
      ! 3/ LOCAL VARIABLES
      ! ------------------
      INTEGER          :: I
      DOUBLE PRECISION :: C1, C2
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
      ! ************************************* !
      ! I - ADIMENSIONAL SKIN STRESS          !
      ! ************************************* !
      c1 = 1.d0/(dens*xmve*grav*dm)
      c2 = 5.d0*sqrt(grav*dens*dm**3)
      CALL os('X=CYZ   ', x=tetap, y=tob, z=mu, c=c1)
!
      CALL os('X=+(Y,C)', x=tetap , y=tetap, c= 1.d-02 )
!
      CALL os('X=Y**C  ', x=ahunzi, y=tetap, c=-1.5d0  )
      CALL os('X=CX    ', x=ahunzi, c= 0.011d0)
      CALL os('X=X+C   ', x=ahunzi, c=-0.3d0  )
!
! BEWARE: AHUNZI CAN BECOME SO LARGE THAT THE HIDING FACTOR BECOMES
! INFINITE; HUNZIKER HIMSELF SUGGESTS THAT IT BE CAPPED TO 2.3.
! THIS INITIALLY ADOPTS A LIMIT OF APPROXIMATELY 10.
! (WHICH IS APPARENT IN TETAP BEING SET TO VALUES .GE. 0.01)
!
!     TODO: REMARK BY JMH: I WOULD STRONGLY RECOMMEND A SINGLE LOOP
!                    WITH THE WHOLE FORMULA, INSTEAD OF PILING
!                    UP CALLS TO OS
!
      DO i = 1, npoin
        hiding%R(i) = (dm/acladm%R(i))**(-ahunzi%R(i))
      ENDDO
      ! ************************************************* !
      ! IV - CORRECTS THE ADIMENSIONAL CRITICAL STRESS    !
      ! ************************************************* !
      CALL os('X=Y/Z   ', x=acp, y=unladm, z=acladm)
      CALL os('X=Y**C  ', x=acp, y=acp   , c=0.33d0)
      CALL os('X=CX    ', x=acp, c=ac)
      ! ********************* !
      ! V - TRANSPORT RATE    !
      ! ********************* !
      CALL os('X=Y-Z   ', x=qsc, y=tetap , z=acp )
      CALL os('X=+(Y,C)', x=qsc, y=qsc   , c=0.d0)
      CALL os('X=XY    ', x=qsc, y=hiding)
      CALL os('X=Y**C  ', x=qsc, y=qsc   , c=1.5d0)

      CALL os('X=CX    ', x=qsc, c=c2)
!======================================================================!
!======================================================================!
      RETURN
      END