bedload_bailard.f

Go to the documentation of this file.

!                   ****************************
                    SUBROUTINE bedload_bailard !
!                   ****************************
!
     &(u2d,v2d,ucmoy,tob,tobw,thetaw,uw,fw,cf,npoin,pi,
     & xmve,grav,dens,xwc,alphaw,qscx,qscy,qssx,qssy,
     & uc3x,uc3y,us4x,us4y,thetac,fcw,qsc,qss,houle)
!
!***********************************************************************
! SISYPHE   V6P1                                   21/07/2011
!***********************************************************************
!
!brief    BAILARD FORMULATION.
!
!history  C. VILLARET (LNHE)
!+        01/10/2003
!+
!+
!
!history  JMH
!+        21/12/2006
!+        V6P0
!+   BEDLOAD_CALCBAIL DELETED; 'HOULE' ARGUMENT ADDED
!
!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 OF WAVES WITH OX
!| CF             |-->| QUADRATIC FRICTION COEFFICIENT
!| DENS           |-->| RELATIVE SENSITY OF SEDIMENT
!| FCW            |---| WAVE-CURRENT FRICTION ANGLE
!| FW             |-->| WAVE FRICTION FACTOR
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| HOULE          |-->| LOGICAL, FOR WAVE EFFECTS
!| NPOIN          |-->| NUMBER OF POINTS
!| PI             |-->| PI
!| QSC            |<->| BEDLOAD TRANSPORT RATE
!| QSCX           |<->| BEDLOAD TRANSPORT RATE IN THE X-DIRECTION
!| QSCY           |<->| BEDLOAD TRANSPORT RATE IN THE Y-DIRECTION
!| QSS            |<->| SUSPENDED LOAD TRANSPORT RATE
!| QSSX           |<->| SUSPENDED LOAD TRANSPORT RATE IN THE X-DIRECTION
!| QSSY           |<->| SUSPENDED LOAD TRANSPORT RATE IN THE Y-DIRECTION
!| THETAC         |<->| CURRENT ANGLE TO THE X DIRECTION
!| THETAW         |-->| ANGLE BETWEEN WAVE AND CURRENT
!| TOB            |<->| BED SHEAR STRESS (TOTAL FRICTION)
!| TOBW           |-->| WAVE INDUCED SHEAR STRESS
!| U2D            |<->| MEAN FLOW VELOCITY X-DIRECTION
!| UC3X           |<->| WORK ARRAY
!| UC3Y           |<->| WORK ARRAY
!| UCMOY          |-->| MEAN CURRENT
!| US4X           |<->| WORK ARRAY
!| US4Y           |<->| WORK ARRAY
!| UW             |-->| ORBITAL WAVE VELOCITY
!| V2D            |<->| MEAN FLOW VELOCITY Y-DIRECTION
!| XMVE           |-->| FLUID DENSITY
!| XWC            |-->| SETTLING VELOCITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,ex_bedload_bailard => bedload_bailard
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
! 2/ GLOBAL VARIABLES
!
      TYPE(bief_obj),   INTENT(IN)    :: U2D,V2D,UCMOY, TOB
      TYPE(bief_obj),   INTENT(IN)    :: TOBW, THETAW, UW, FW, CF
      INTEGER,          INTENT(IN)    :: NPOIN
      LOGICAL,          INTENT(IN)    :: HOULE
      DOUBLE PRECISION, INTENT(IN)    :: PI, XMVE, GRAV, DENS, XWC
      TYPE(bief_obj),   INTENT(INOUT) :: ALPHAW        ! WORK ARRAY T1
      TYPE(bief_obj),   INTENT(INOUT) :: QSCX, QSCY    ! WORK ARRAY T2 AND T3
      TYPE(bief_obj),   INTENT(INOUT) :: QSSX, QSSY    ! WORK ARRAY T4 AND T5
      TYPE(bief_obj),   INTENT(INOUT) :: UC3X, UC3Y    ! WORK ARRAY T6 AND T7
      TYPE(bief_obj),   INTENT(INOUT) :: US4X, US4Y    ! WORK ARRAY T8 AND T9
      TYPE(bief_obj),   INTENT(INOUT) :: THETAC, FCW   ! WORK ARRAY T10 AND T11
      TYPE(bief_obj),   INTENT(INOUT) :: QSC, QSS
!
! 3/ LOCAL VARIABLES
!
      INTEGER                     :: I
      DOUBLE PRECISION            :: C3, C4, PHI
      DOUBLE PRECISION, PARAMETER :: EPSC = 0.21d0   ! BEDLOAD
      DOUBLE PRECISION, PARAMETER :: EPSS = 0.025d0  ! SUSPENSION
      DOUBLE PRECISION            :: U3X, U3Y, NUM
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
!     CASE WITH WAVES
!
      IF(houle) THEN
!
!     ANGLE OF VELOCITY WITH OX (IN RADIANS)
!
      CALL bedload_direction(u2d,v2d,npoin,pi,thetac)
!
!     ANGLE OF WAVES WITH OX (IN RADIANS)
!
      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)
!
!     PARAMETERS  ,
!
!
!     US4X AND US4Y ARE WORK ARRAYS, THEIR STRUCTURE IS GIVEN HERE
!     THE STRUCTURE OF THETAC (CATHERINE DON'T REMOVE THIS PLEASE)
      CALL cpstvc(thetac,us4x)
      CALL cpstvc(thetac,us4y)
!
      DO i = 1, npoin
        ! ********************* !
        ! I - CURRENT REFERENCE SYSTEM !
        ! ********************* !
        u3x = ucmoy%R(i)**3
     &      + ucmoy%R(i)*uw%R(i)**2 * (1 + cos(2.d0*alphaw%R(i))/ 2.d0)
        u3y = ucmoy%R(i)*uw%R(i)**2 * sin(2.d0*alphaw%R(i)) / 2.d0
        ! ********************************************** !
        ! II - 3RD ORDER MOMENTUM (LINEAR WAVE THEORY)   !
        ! ********************************************** !
        uc3x%R(i) = u3x * cos(thetac%R(i)) - u3y * sin(thetac%R(i))
        uc3y%R(i) = u3x * sin(thetac%R(i)) + u3y * cos(thetac%R(i))
        ! ************************************************************ !
        ! III -  4TH ORDER MOMENTUM (COLINEAR WAVES AND CURRENTS)      !
        ! ************************************************************ !
        num = ( 8.d0*ucmoy%R(i)**4 + 3.d0*uw%R(i)**4
     &          + 24.d0*(ucmoy%R(i)**2)*(uw%R(i)**2) )*0.125d0
        us4x%R(i) = num * cos(thetac%R(i))
        us4y%R(i) = num * sin(thetac%R(i))
      ENDDO
      ! *********************************************** !
      ! IV -  FRICTION COEFFICIENT WAVE + CURRENT       !
      ! *********************************************** !
      CALL bedload_interact
     &     (ucmoy,tobw,tob,alphaw,fw,cf,uw,npoin,xmve,fcw)
      ! ******************************** !
      ! V - TRANSPORT RATES              !
      ! ******************************** !
      phi = pi   / 6.d0  ! FRICTION ANGLE
      c3  = epsc / (grav*dens*tan(phi))
      c4  = epss / (grav*dens*xwc)
      CALL os('X=CYZ   ', x=qscx, y=fcw,  z=uc3x, c=c3)
      CALL os('X=CYZ   ', x=qscy, y=fcw,  z=uc3y, c=c3)
      CALL os('X=CYZ   ', x=qssx, y=fcw,  z=us4x, c=c4)
      CALL os('X=CYZ   ', x=qssy, y=fcw,  z=us4y, c=c4)
!
!     CASE WITHOUT WAVES
!
      ELSE
!
        WRITE(lu,*) 'BAILARD WITHOUT WAVES NOT PROGRAMMED'
        CALL plante(1)
        stop
!
      ENDIF
!
!     NORMS
!
      CALL os('X=N(Y,Z)', x=qsc,  y=qscx, z=qscy)
      CALL os('X=N(Y,Z)', x=qss,  y=qssx, z=qssy)
!======================================================================!
!======================================================================!
      RETURN
      END