bedload_effpnt.f

Go to the documentation of this file.

!                   *************************
                    SUBROUTINE bedload_effpnt
!                   *************************
!
     & (maskel,liqbor,s,zf,npoin,nptfr,ielmt,kent,
     &  beta,pi,msk,mesh,dzfdx,dzfdy,cteta,steta,
     &  coef,calfa,salfa,slopeff,phised,devia,beta2,
     &  tob,xmvs,xmve,dm,grav,unsv2d)
!
!***********************************************************************
! SISYPHE   V6P3                                   12/02/2013
!***********************************************************************
!
!brief    COMPUTES THE PARAMETERS OF THE SLOPE EFFECT.
!
!history  E. PELTIER; C. LENORMANT; J.-M. HERVOUET
!+        11/09/1995
!+        V5P1
!+
!
!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  Pablo Tassi PAT (EDF-LNHE)
!+        12/02/2013
!+        V6P3
!+ Correction by Rebekka Kopmann (BAW):
!+ Avoiding the slope effect at the open boundaries
!+ So far only the magnitude is set to 1, but the directions are due to the slope effect.
!+ Specially at the outlet boundary, this can create problems.
!+ If a bar is moved out of the model there can be the situation,
!+ that you have a negative slope in longitudinal, which creates in inflow of the bed load
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| BETA           |-->| COEFFICIENT FOR SLOPING BED EFFECT ( KOCH AND FLOKSTRA)
!| BETA2          |-->| COEFFICIENT FOR THE DEVIATION  (TALMON ET AL.)
!| CALFA          |<->| COSINUS OF THE ANGLE BETWEEN TRANSPORT RATE AND X-AXIS
!| COEF           |<->| CORRECTION COEFFICIENT FOR THE INTENSITY OF BEDLOAD TRANSPORT
!| CTETA          |<->| COSINUS OF THE ANGLE BETWEEN MEAN FLOW AND X-AXIS
!| DEVIA          |-->| SLOPE EFFECT FORMULA FOR DEVIATION
!| DM             |-->| SEDIMENT GRAIN DIAMETER
!| DZFDX          |<->| BOTTOM SLOPE IN THE X-DIRECTION
!| DZFDY          |<->| BOTTOM SLOPE IN THE Y-DIRECTION
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| IELMT          |-->| NUMBER OF ELEMENTS
!| KENT           |-->| CONVENTION FOR LIQUID INPUT WITH PRESCRIBED VALUE
!| LIQBOR         |-->| TYPE OF BOUNDARY CONDITION FOR QS
!| MASKEL         |-->| MASKING OF ELEMENTS
!| MESH           |<->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| PHISED         |-->| ANGLE OF REPOSE OF THE SEDIMENT
!| PI             |-->| PI
!| S              |-->| VOID STRUCTURE
!| SALFA          |<->| SINUS OF THE ANGLE BETWEEN TRANSPORT RATE AND X-AXIS
!| SLOPEFF        |-->| LOGICAL, SLOPING BED EFFECT OR NOT
!| STETA          |<->| COSINUS OF THE ANGLE BETWEEN MEAN FLOW AND Y-AXIS
!| TOB            |<->| BED SHEAR STRESS (TOTAL FRICTION)
!| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
!| XMVE           |-->| FLUID DENSITY
!| XMVS           |-->| WATER DENSITY
!| ZF             |-->| ELEVATION OF BOTTOM
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,ex_bedload_effpnt => bedload_effpnt
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
      ! 2/ GLOBAL VARIABLES
      ! -------------------
      TYPE(bief_obj),   INTENT(IN)    :: MASKEL,LIQBOR,S,UNSV2D
      TYPE(bief_obj),   INTENT(IN)    :: ZF, TOB
      INTEGER,          INTENT(IN)    :: NPOIN, NPTFR, IELMT, KENT
      INTEGER,          INTENT(IN)    :: SLOPEFF,DEVIA
      DOUBLE PRECISION, INTENT(IN)    :: BETA, PI, PHISED, BETA2
      DOUBLE PRECISION, INTENT(IN)    :: XMVS, XMVE, GRAV, DM
      LOGICAL,          INTENT(IN)    :: MSK
      TYPE(bief_mesh),  INTENT(INOUT) :: MESH
      TYPE(bief_obj),   INTENT(INOUT) :: DZFDX, DZFDY
      TYPE(bief_obj),   INTENT(INOUT) :: CTETA,STETA
      TYPE(bief_obj),   INTENT(INOUT) :: COEF, CALFA, SALFA
!
      ! 3/ LOCAL VARIABLES
      ! ------------------
      INTEGER          :: I, K
      DOUBLE PRECISION :: C,ZETA,C1,CALPHA,SALPHA,AA,BB
      DOUBLE PRECISION :: CPSI,SPSI,DZF,TANPHI,CZETA,SZETA,SURBETA2
      DOUBLE PRECISION :: NORM,TT1
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
!
!     TODO : IF WE SWAP THE DEVIA AND THE SLOPEFF ACTIONS BELOW
!            IT IS PROBABLY NOT USEFUL TO DO A COPY OF CALFA AND SALFA ?
!
!     COS AND SIN TETA (CALFA AND SALFA HAVE ALREADY BEEN COMPUTED
!                       IN BEDLOAD_SOLIDISCHARGE, SO JUST A COPY HERE)
!     TETA = ANGLE OF THE FLOW WITH OX
!
      CALL os('X=Y     ',x=cteta,y=calfa)
      CALL os('X=Y     ',x=steta,y=salfa)
!
!----------------------------------------------------------------------
!
!     COMPUTES THE SLOPE  : D(ZF)/DX ET D(ZF)/DY (AT THE NODES)
!
      CALL vector(dzfdx, '=', 'GRADF          X',ielmt,1.d0,zf,s,s,
     &            s,s,s,mesh,msk,maskel)
      CALL vector(dzfdy, '=', 'GRADF          Y',ielmt,1.d0,zf,s,s,
     &            s,s,s,mesh,msk,maskel)
!
      IF(ncsize.GT.1) THEN
        CALL parcom(dzfdx,2,mesh)
        CALL parcom(dzfdy,2,mesh)
      ENDIF
!
      CALL os('X=XY    ',x=dzfdx,y=unsv2d)
      CALL os('X=XY    ',x=dzfdy,y=unsv2d)
!
!======================================================================
!
!     COMPUTES THE SOLID TRANSPORT ANGLE ALFA = TETA + DEVIATION
!
!     1 : KOCH AND FLOKSTRA
!
      IF(devia==1) THEN
!
      c = 2.d0*(xmvs-xmve)*grav*dm/3.d0
      DO i=1,npoin
        tt1=c/max(tob%R(i),1.d-10)
        aa=steta%R(i)-tt1*dzfdy%R(i)
        bb=cteta%R(i)-tt1*dzfdx%R(i)
        norm=max(sqrt(aa**2+bb**2),1.d-10)
        salfa%R(i)=aa/norm
        calfa%R(i)=bb/norm
      ENDDO
!
!     2 : TALMON ET AL. JHR 1995 33(4)
!
      ELSEIF(devia==2) THEN
!
      surbeta2=1.d0/beta2
      c = (xmvs-xmve)*grav*dm*surbeta2**2
      DO i=1,npoin
        tt1=sqrt(c/max(tob%R(i),1.d-10))
        aa=steta%R(i)-tt1*dzfdy%R(i)
        bb=cteta%R(i)-tt1*dzfdx%R(i)
        norm=max(sqrt(aa**2+bb**2),1.d-10)
        salfa%R(i)=aa/norm
        calfa%R(i)=bb/norm
      ENDDO
!
      ENDIF
!
!======================================================================
!
!     COMPUTES THE COEFFICIENT TO TAKE THE SLOPE EFFECT ON THE MAGNITUDE
!     OF THE SOLID TRANSPORT INTO ACCOUNT
!
!     METHOD 1 (EMPIRICAL)
!
      IF(slopeff==1) THEN
!
        DO i=1,npoin
          coef%R(i)=max(0.d0,
     &    1.d0-beta*(dzfdx%R(i)*cteta%R(i)+dzfdy%R(i)*steta%R(i)) )
        ENDDO
!
!     METHOD 2 : SOULSBY 1997 DYNAMICS OF MARINE SANDS P107-108
!
      ELSEIF(slopeff.EQ.2) THEN
!
        tanphi = tan(phised*pi/180.d0)
!
        DO i=1,npoin
!
!         COSINE AND SINE OF THE DIRECTION OF THE SLOPE
          dzf=sqrt(dzfdx%R(i)**2+dzfdy%R(i)**2)
          IF(dzf.GT.1.d-12) THEN
            calpha=dzfdx%R(i)/dzf
            salpha=dzfdy%R(i)/dzf
          ELSE
            calpha=1.d0
            salpha=0.d0
          ENDIF
!
!         ZETA: ANGLE OF THE SLOPE WITH HORIZONTAL (BETA IN SOULSBY)
          zeta=atan(dzf)
          czeta=cos(zeta)
          szeta=sin(zeta)
!
!         PSI: ANGLE OF THE CURRENT WITH THE SLOPE DIRECTION
!         PSI=TETA%R(I)-ALPHA
          cpsi=cteta%R(i)*calpha+steta%R(i)*salpha
          spsi=steta%R(i)*calpha-cteta%R(i)*salpha
          c1=(czeta*tanphi)**2-(spsi*szeta)**2
          coef%R(i)=max((cpsi*szeta+sqrt(max(c1,0.d0)))/tanphi,0.d0)
          coef%R(i)=max(coef%R(i),0.d0)
!
        ENDDO
!
      ENDIF
!
! ********************************************************************* !
!     V - TREATS THE BOUNDARY POINTS WITH IMPOSED RATE                  !
!         QS IS NOT MODIFIED WHEN SPECIFIED BY THE USER                 !
! ********************************************************************* !
!
! RK no slope effect at the open boundaries
      DO k = 1, nptfr
        IF (liqbor%I(k) == kent) THEN
          coef%R(mesh%NBOR%I(k)) = 1.d0
          calfa%R(mesh%NBOR%I(k)) = cteta%R(mesh%NBOR%I(k))
          salfa%R(mesh%NBOR%I(k)) = steta%R(mesh%NBOR%I(k))
        ENDIF
!     R.K. MAY 2007
!     KSORT = 4
        IF (liqbor%I(k) == 4) THEN
          coef%R(mesh%NBOR%I(k)) = 1.d0
          calfa%R(mesh%NBOR%I(k)) = cteta%R(mesh%NBOR%I(k))
          salfa%R(mesh%NBOR%I(k)) = steta%R(mesh%NBOR%I(k))
        ENDIF
      ENDDO
!
!======================================================================
!======================================================================
!
      RETURN
      END SUBROUTINE bedload_effpnt