bedload_dibwat.f

Go to the documentation of this file.

!                   ***************************
                    SUBROUTINE bedload_dibwat !
!                   ***************************
!
     &(u2d,v2d,ucmoy, cf, tob, tobw, uw, tw, fw, thetaw, npoin,
     & xmve, dens, grav, dm, xwc, pi, alphaw, t2, t3, ucw, ucn,
     & uw1, uw2, tw1, tw2, thetac, fcw, qsc)
!
!***********************************************************************
! SISYPHE   V6P1                                   21/07/2011
!***********************************************************************
!
!brief    DIBAJNIA & WATANABE FORMULATION (1992).
!
!history  C. VILLARET (LNHE)
!+        15/11/2003
!+
!+
!
!history  JMH
!+        09/05/2007
!+        V5P7
!+   CHECKS FOR GAMMA=0 DIVISION
!
!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
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ALPHAW         |<->| ANGLE BETWEEN WAVE AND CURRENT
!| CF             |-->| QUADRATIC FRICTION COEFFICIENT
!| DENS           |-->| RELATIVE DENSITY
!| DM             |-->| SEDIMENT GRAIN DIAMETER
!| FCW            |-->| WAVE-CURRENT FRICTION FACTOR
!| FW             |-->| WAVE FRICTION FACTOR
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| NPOIN          |-->| NUMBER OF POINTS
!| PI             |-->| PI
!| QSC            |<->| BED LOAD TRANSPORT
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T3             |<->| WORK BIEF_OBJ STRUCTURE
!| THETAC         |<->| CURRENT ANGLE TO THE X AXIS
!| THETAW         |-->| ANGLE BETWEEN WAVE AND CURRENT
!| TOB            |-->| BED SHEAR STRESS (TOTAL FRICTION)
!| TOBW           |-->| WAVE INDUCED SHEAR STRESS
!| TW             |-->| WAVE PERIOD
!| TW1            |<->| MID PERIOD (U(T)>0)
!| TW2            |<->| MID PERIOD (U(T)<0)
!| U2D            |<->| MEAN FLOW VELOCITY X-DIRECTION
!| UCMOY          |-->| MEAN CURRENT
!| UCN            |<->| MEAN CURRENT AT AN ANGLE TO THE WAVE
!| UCW            |<->| MEAN CURRENT PROJECTED IN THE WAVE DIRECTION
!| UW             |-->| ORBITAL WAVE VELOCITY
!| UW1            |<->| MEAN CURRENT IN THE WAVE DIRECTION
!| UW2            |<->| MEAN CURRENT IN THE OPPOSITE DIRECTION
!| V2D            |<->| MEAN FLOW VELOCITY Y-DIRECTION
!| XMVE           |-->| FLUID DENSITY
!| XWC            |-->| SETTLING VELOCITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,ex_bedload_dibwat => bedload_dibwat
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
      ! 2/ GLOBAL VARIABLES
      ! -------------------
      TYPE(bief_obj),   INTENT(IN)    :: U2D,V2D,UCMOY, CF, TOB, TOBW
      TYPE(bief_obj),   INTENT(IN)    :: UW, TW, FW
      TYPE(bief_obj),   INTENT(IN)    :: THETAW
      INTEGER,          INTENT(IN)    :: NPOIN
      DOUBLE PRECISION, INTENT(IN)    :: XMVE, DENS, GRAV, DM, XWC, PI
      TYPE(bief_obj),   INTENT(INOUT) :: ALPHAW          ! WORK ARRAY T1
      TYPE(bief_obj),   INTENT(INOUT) :: T2, T3          !
      TYPE(bief_obj),   INTENT(INOUT) ::  UCW, UCN ! WORK ARRAY T4, T5, T6
      TYPE(bief_obj),   INTENT(INOUT) :: UW1, UW2, TW1   ! WORK ARRAY T7, T8, T9
      TYPE(bief_obj),   INTENT(INOUT) :: TW2, THETAC     ! WORK ARRAY T10, T11
      TYPE(bief_obj),   INTENT(INOUT) :: FCW, QSC        ! WORK ARRAY T12
!
      ! 3/ LOCAL VARIABLES
      ! ------------------
      INTEGER                     :: I
      DOUBLE PRECISION            :: UCW2, SHIELDS, OMEGACR, OMEGA1
      DOUBLE PRECISION            :: OMEGA2, W1, WP1, W2, WP2, GAMMAW
      DOUBLE PRECISION            :: GAMMAN, GAMMA, QSW, QSN
      DOUBLE PRECISION, PARAMETER :: ZERO  = 1.d-6
      DOUBLE PRECISION, PARAMETER :: ALPHA = 0.0001d0
      DOUBLE PRECISION, PARAMETER :: BETA  = 0.55d0
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
!     ANGLE OF VELOCITY WITH OX
!
      CALL bedload_direction(u2d,v2d,npoin,pi,thetac)
!
!     PROJECTIONS OF UCMOY IN THE WAVE DIRECTION
!
      CALL os('X=CY    ', x=alphaw, y=thetaw, c=-pi/180.d0)
      CALL os('X=X+C   ', x=alphaw, c=0.5d0*pi)
      CALL os('X=Y-Z   ', x=alphaw, y=alphaw, z=thetac)
      CALL os('X=COS(Y)', x=t2    , y=alphaw)
      CALL os('X=SIN(Y)', x=t3    , y=alphaw)
      CALL os('X=YZ    ', x=ucw   , y=ucmoy , z=t2)
      CALL os('X=CYZ   ', x=ucn   , y=ucmoy , z=t3, c= -1.d0)
!
!     QUADRATIC VELOCITIES (UW1=U(T)> 0 AND UW2=U(T)
!
      CALL bedload_calcdw(ucw,uw,tw,npoin,pi,uw1,uw2,tw1,tw2)
!
!     FRICTION COEFFICIENT
!
      CALL bedload_interact
     &     (ucmoy,tobw,tob,alphaw,fw,cf,uw,npoin,xmve,fcw)
!
      DO i = 1, npoin
!
!       SHIELDS PARAMETER
        ucw2    = ucmoy%R(i)**2 + 0.5d0 * uw%R(i)**2
        shields = fcw%R(i)*ucw2/dens/grav/dm
!
!       CRITICAL OMEGA (TEMPERVILLE AND GUIZA,2000)  !
        IF (shields <= 0.2d0) THEN
          omegacr = 0.03d0
        ELSEIF (shields <= 0.4d0) THEN
          omegacr = 1.d0 - sqrt(1.d0-((shields-0.2d0)/0.58d0)**2)
        ELSEIF (shields <= 1.5d0) THEN
          omegacr = 0.8104d0 * sqrt(shields) - 0.4225d0
        ELSE
          omegacr = 0.7236d0 * sqrt(shields) - 0.3162d0
        ENDIF
!
!       OMEGA1, OMEGA2  !
        IF(tw1%R(i) > zero) THEN
          omega1 = uw1%R(i)**2 / (2.d0*dens*grav*xwc*tw1%R(i))
        ELSE
          omega1 = zero
        ENDIF
!
        IF(tw2%R(i) > zero) THEN
          omega2 = uw2%R(i)**2 / (2.d0*dens*grav*xwc*tw2%R(i))
        ELSE
          omega2 = zero
        ENDIF
!
!       QUANTITIES OF SAND DEPOSITED AND
!       IN SUSPENSION AT EACH PHASE OF THE CYCLE
!
!       PHASE 1
        IF (omega1 <= omegacr) THEN
          w1  = 2.d0 * omega1 * xwc * tw1%R(i) / dm
          wp1 = 0.d0
        ELSE
          w1  = 2.d0 * omegacr * xwc * tw1%R(i) / dm
          wp1 = 2.d0 * (omega1-omegacr) * xwc * tw1%R(i) / dm
        ENDIF
!
!       PHASE 2
        IF (omega2 <= omegacr) THEN
          w2  = 2.d0 * omega2 * xwc * tw2%R(i) / dm
          wp2 = 0.d0
        ELSE
          w2  = 2.d0 * omegacr * xwc * tw2%R(i) / dm
          wp2 = 2.d0 * (omega2-omegacr) * xwc * tw2%R(i) / dm
        ENDIF
!
!       GAMMAW, GAMMAN,GAMMA
        IF ((uw2%R(i)*tw2%R(i) + uw1%R(i)*tw1%R(i)) > 0.d0 ) THEN
          gammaw = (  uw1%R(i) * tw1%R(i) * (w1**3+wp2**3)
     &              - uw2%R(i) * tw2%R(i) * (w2**3+wp1**3))
     &           / ( uw1%R(i)*tw1%R(i)+uw2%R(i)*tw2%R(i))
        ELSE
          gammaw = (2.d0 * ucw%R(i)**2 / dens / grav / dm)**3
        ENDIF
        gamman = (2.d0 * ucn%R(i)**2 / dens / grav / dm)**3
        gamma  = max(sqrt(gammaw**2 + gamman**2),1.d-10)
!
!       SOLID TRANSPORT IN THE WAVE DIRECTION : QSW
!       AND IN THE NORMAL DIRECTION : QSN
!
        qsw      = alpha * gammaw * xwc * dm / gamma**(1.d0-beta)
        qsn      = alpha * gamman * xwc * dm / gamma**(1.d0-beta)
        qsc%R(i) = sqrt(qsw**2 + qsn**2)
!
      ENDDO
!
!----------------------------------------------------------------------
!
      RETURN
      END