bedload_formula.f

Go to the documentation of this file.

!                   **************************
                    SUBROUTINE bedload_formula
!                   **************************
!
     &(u2d,v2d,ucmoy,hn,cf,mu,tob,tobw,uw,tw,thetaw,fw,
     & acladm, unladm,ksp,ksr,ava,npoin,icf,hidfac,xmvs,xmve,
     & dm,grav,vce,hmin,xwc,d90,karman,zero,
     & pi,susp, ac, hiding, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10,
     & t11,tetap, qsc, qss,ielmt,seccurrent,slopeff,
     & coefpn,calfa,salfa,bijk,houle, sanfra)
!
!***********************************************************************
! SISYPHE   V8P0                                   12/09/2018
!***********************************************************************
!
!brief    COMPUTES THE BED-LOAD TRANSPORT.
!
!history  BUI MINH DUC
!+        **/01/2002
!+        V5P2
!+
!
!history  C. VILLARET
!+        **/10/2003
!+        V5P4
!+
!
!history  F. HUVELIN
!+        12/01/2005
!+        V5P6
!+
!
!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
!+
!history  J-M HERVOUET (EDF-LNHE)
!+        27/02/2012
!+        V6P2
!+  ALPHA suppressed, was no longer used
!
!history  R KOPMANN (BAW)
!+        10/05/2016
!+        V7P2
!+ CALFA,SALFA dependent of grain classes
!
!history  F.CORDIER & P.TASSI (EDF-LNHE)
!+        12/09/2018
!+        V8P0
!+  PARSING of SANFRA for the formula of WILCOCK and CROWE (2003)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AC             |<->| SHIELDS PARAMETER
!| ACLADM         |-->| MEAN DIAMETER
!| AVA            |-->| PERCENT AVAILABLE
!| BIJK           |-->| EMPIRICAL COEFFICIENT
!| CF             |-->| QUADRATIC FRICTION COEFFICIENT
!| COEFPN         |-->| COEFFICIENT FOR SLOPING BED EFFECTS
!| D90            |-->| D90
!| DM             |-->| DIAMETER OF THE CLASS
!| FW             |-->| WAVE FRICTION COEFFICIENT
!| GRAV           |-->| GRAVITY
!| HIDFAC         |-->| HIDING FACTOR FORMULA
!| HIDING         |<->| HIDING FACTOR
!| HMIN           |-->| MININMUM WATER DEPTH
!| HN             |-->| WATER DEPTH
!| HOULE          |-->| EFFECT OF WAVE
!| ICF            |-->| CHOICE OF FORMULA
!| IELMT          |-->| NUMBER OF ELEMENTS
!| KARMAN         |-->| VON KARMAN COEFFICIENT
!| KSP            |-->| SKIN BED ROUGHNESS (M)
!| KSR            |-->| RIPPLE BED ROUGHNESS (M)
!| MU             |-->| CORRECTION FOR SKIN FRICTION
!| NPOIN          |-->| NUMBER OF POINTS
!| PI             |-->| PI
!| QSC            |<->| BED LOAD TRANSPORT RATE (m2/S)
!| QSS            |<->| SUSPENDED LOAD TRANSPORT RATE (M2/S)
!| SANFRA         |-->| SAND FRACTION (WILCOCK AND CROWE, 2003)
!| SECCURRENT     |-->| EFFECT OF SECUNDARY CURRENTS
!| SLOPEFF        |-->| FORMULA FOR SLOPING BED EFFECTS
!| SUSP           |-->| SUSPENSION TREATMENT
!| T1             |<->| WORKING ARRAYS
!| T10            |<->| --
!| T11            |<->| --
!| T2             |<->| --
!| T3             |<->| --
!| T4             |<->| --
!| T5             |<->| --
!| T6             |<->| --
!| T7             |<->| --
!| T8             |<->| --
!| T9             |<->| --
!| TETAP          |<->| ADIMENSIONAL SKIN FRICTION
!| THETAW         |-->| WAVE/CURRENT ANGLE
!| TOB            |-->| TOTAL BED SHEAR STRESS (N/M2)
!| TOBW           |-->| WAVE INDUCED BED SHEAR STRESS (N/M2)
!| TW             |-->| WAVE PERIOD (S)
!| U2D            |-->| LONGITUDINAL VELOCITY (m/S)
!| UCMOY          |-->| CURRENT INTENSITY (M/S)
!| UNLADM         |-->| DIAMETER OF LAYER 2 (M)
!| UW             |-->| WAVE ORBITAL VELOCITY (M/S)
!| V2D            |-->| TRANSVERSAL VELOCITY (M/S)
!| VCE            |-->| FLUID KINEMATIC VISCOSITY (M2/S)
!| XMVE           |-->| FLUID DENSITY (KG/M3)
!| XMVS           |-->| SEDIMENT DENSITY (KG/M3)
!| XWC            |-->| SETTLING VELOCITY (M/S)
!| ZERO           |-->| ZERO
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,ex_bedload_formula => bedload_formula
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj),   INTENT(IN)    :: U2D, V2D, UCMOY,HN, CF, TOB
      TYPE(bief_obj),   INTENT(IN)    :: MU,TOBW, UW, TW, THETAW, FW
      TYPE(bief_obj),   INTENT(IN)    :: ACLADM,UNLADM,KSR,KSP
      INTEGER,          INTENT(IN)    :: NPOIN, ICF, HIDFAC,IELMT
      DOUBLE PRECISION, INTENT(IN)    :: XMVS, XMVE, DM, GRAV, VCE
      DOUBLE PRECISION, INTENT(IN)    :: HMIN, XWC, D90
      DOUBLE PRECISION, INTENT(IN)    :: KARMAN, ZERO, PI
      LOGICAL,          INTENT(IN)    :: SUSP,SECCURRENT,HOULE
      DOUBLE PRECISION, INTENT(INOUT) :: AC
      TYPE(bief_obj),   INTENT(INOUT) :: HIDING
      TYPE(bief_obj),   INTENT(INOUT) :: T1, T2, T3, T4, T5, T6, T7
      TYPE(bief_obj),   INTENT(INOUT) :: T8, T9, T10,T11
      TYPE(bief_obj),   INTENT(INOUT) :: TETAP ! WORK ARRAY T12
      TYPE(bief_obj),   INTENT(INOUT) :: QSC, QSS
      TYPE(bief_obj),   INTENT(INOUT) :: COEFPN, CALFA, SALFA
      INTEGER,          INTENT(IN)    :: SLOPEFF
!
      DOUBLE PRECISION, INTENT (IN) :: BIJK,AVA(npoin)
      DOUBLE PRECISION, INTENT (IN) :: SANFRA(npoin)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER                     :: I
      DOUBLE PRECISION            :: DENS,DSTAR
      DOUBLE PRECISION, PARAMETER :: ZERO_LOCAL = 1.d-6
      DOUBLE PRECISION            :: C1

!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
      ! *************************** !
      ! I - ADIMENSIONAL PARAMETERS !
      ! *************************** !
      dens  = (xmvs - xmve )/ xmve
      dstar = dm*(grav*dens/vce**2)**(1.d0/3.d0)
      ! ************************ !
      ! II -  SKIN FRICTION      !
      ! ************************ !
!
      c1 = 1.d0/(dens*xmve*grav*dm)
      CALL os('X=CYZ   ', x=tetap, y=tob,z=mu,  c=c1)
      CALL os('X=+(Y,C)', x= tetap,y=tetap, c=zero_local)
!
      IF(seccurrent) CALL bedload_seccurrent(ielmt,calfa,salfa)
      ! ****************************************** !
      ! IV - COMPUTES 2 TRANSPORT TERMS            !
      !      QSS : SUSPENSION                      !
      !      QSC : BEDLOAD                         !
      ! ****************************************** !
      ! ===================================== !
      ! IV(1) - MEYER-PETER-MULLER FORMULATION!
      !         FOR BEDLOAD ONLY              !
      ! ===================================== !
      IF(icf == 1) THEN
!
        CALL bedload_meyer(tetap,hiding,hidfac,dens,grav,dm,ac,
     &                     t1,qsc,slopeff,coefpn)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)
        ENDDO
!
      ! =========================== !
      ! IV(2) - EINSTEIN FORMULATION!
      !         FOR BEDLOAD ONLY    !
      ! =========================== !
      ELSEIF(icf == 2) THEN
        CALL bedload_einst(tetap,npoin,dens,grav,dm,dstar,qsc)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! =================================== !
      ! IV(30) - ENGELUND-HANSEN FORMULATION!
      !          FOR TOTAL TRANSPORT        !
      ! =================================== !
      ELSEIF(icf == 30) THEN
!       V6P0 MU IS USED INSTEAD OF CF
!       BEWARE: DIFFERENCES
!       CALL BEDLOAD_ENGEL(TETAP,DENS,GRAV,DM,QSC)
!       BACK TO EARLIER VERSION OF BEDLOAD_ENGEL
        CALL bedload_engel(tob,cf,dens,grav,dm,xmve,qsc)
!       ARBITRARY DISTRIBUTION
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! ======================================== !
      ! IV(3) - ENGELUND-HANSEN FORMULATION      !
      !         MODIFIED: CHOLLET ET CUNGE       !
      !         FOR TOTAL TRANSPORT              !
      ! ======================================== !
      ELSEIF(icf == 3) THEN
!       KSP IS USED INSTEAD OF CFP
        CALL bedload_engel_cc
     &       (tetap,cf,npoin,grav,dm,dens,t1,qsc)
!       ARBITRARY DISTRIBUTION
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! ============================== !
      ! IV(4) - BIJKER FORMULATION     !
      !         FOR BEDLOAD + SUSPENSION !
      ! ============================== !
      ELSEIF (icf == 4) THEN
        CALL bedload_bijker
     &   (tobw,tob,mu,ksp,ksr,hn,npoin,dm,dens,xmve,grav,
     &    xwc,karman,zero,t4,t7,t8,t9,qsc,qss,bijk,houle)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
          qss%R(i)=qss%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! ============================== !
      ! IV(5) - SOULSBY FORMULATION    !
      !         FOR BEDLOAD + SUSPENSION !
      ! ============================== !
      ELSEIF (icf == 5) THEN
        CALL bedload_soulsby
     &       (ucmoy,hn,uw,npoin,dens,grav,dm,dstar,
     &        d90,qsc,qss)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
          qss%R(i)=qss%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! ================================================== !
      ! IV(6) - HUNZIKER / MEYER-PETER & MULLER FORMULATION!
      !         FOR BEDLOAD ONLY                           !
      ! ================================================== !
      ELSEIF (icf == 6) THEN
        CALL bedload_hunz_meyer
     &       (tob, mu, acladm, unladm, npoin, dens, xmve, grav,
     &        dm, ac, t1, t2, t3, hiding, qsc)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)
        ENDDO
!
      ! =========================== !
      ! IV(7) - VAN RIJN FORMULATION!
      !         FOR BEDLOAD ONLY    !
      ! =========================== !
      ELSEIF (icf == 7) THEN
!
        CALL bedload_vanrijn
     &       (tetap,npoin,dm,dens,grav,dstar,ac,qsc)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! ============================== !
      ! IV(8) - BAILARD FORMULATION    !
      !         FOR BEDLOAD + SUSPENSION !
      ! ============================== !
      ELSEIF (icf == 8) THEN
!
        CALL bedload_bailard
     &       (u2d,v2d,ucmoy,tob,tobw,thetaw,uw,fw,cf,npoin,
     &        pi,xmve,grav,dens,xwc,t1,t2,t3,t4,t5,t6,t7,
     &        t8,t9,t10,t11,qsc,qss,houle)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
          qss%R(i)=qss%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! ======================================= !
      ! IV(9) - DIBAJNIA AND WATANABE FORMULATION!
      !         FOR TOTAL TRANSPORT             !
      ! ======================================= !
      ELSEIF(icf == 9) THEN
!
        CALL bedload_dibwat
     &       (u2d,v2d,ucmoy, cf, tob, tobw, uw, tw, fw, thetaw,
     &        npoin, xmve, dens, grav, dm, xwc, pi, t1, t2, t3, t4,
     &        t5, t6, t7, t8, t9, t10, t11, qsc)
!       ARBITRARY DISTRIBUTION
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
        ENDDO
!
      ! ======================================= !
      ! IV(10) - WILCOCK AND CROWE 2003 FORMULATION!
      !         NON-UNIFORM TRANSPORT           !
      ! ======================================= !
      ELSEIF(icf == 10) THEN
!
        CALL bedload_wilcock_crowe
     & (tob, mu, acladm, dm, ava, grav, xmve, xmvs, sanfra, qsc, ac, t1,
     &  slopeff, coefpn)
!
      ! ============================================ !
      ! IV(0) - USER-DEFINED FORMULATION             !
      ! ============================================ !
      ELSEIF (icf == 0) THEN
        CALL qsform
     &       (u2d, v2d, tob, hn, xmve, tetap, mu, npoin, dm,
     &       dens, grav, dstar, ac, qsc, qss)
        DO i=1,npoin
          qsc%R(i)=qsc%R(i)*ava(i)*hiding%R(i)
          qss%R(i)=qss%R(i)*ava(i)*hiding%R(i)
        ENDDO
      ! ================= !
      ! IV(ELSE) - ERROR  !
      ! ================= !
      ELSE
        WRITE(lu,201) icf
201     FORMAT(1x,'TRANSP : TRANSPORT FORMULA UNKNOWN:',1i6)
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     WHEN SUSPENSION IS NOT ASKED SPECIFICALLY, SOME BEDLOAD TRANSPORT
!     FORMULAS GIVE A VALUE
!
      IF(.NOT.susp) THEN
        IF(icf.EQ.4.OR.icf.EQ.5.OR.icf.EQ.8.OR.icf.EQ.0) THEN
          DO i = 1,npoin
            qsc%R(i) = qsc%R(i) + qss%R(i)
          ENDDO
        ELSE
!         TODO: NOTE JMH: IS THIS REALLY USEFUL ???
          DO i = 1,npoin
            qss%R(i) = 0.d0
          ENDDO
        ENDIF
      ENDIF
!
!=======================================================================
!=======================================================================
!
      RETURN
      END