bedload_nerbed_vf.f

Go to the documentation of this file.

!                   ******************************
                    SUBROUTINE bedload_nerbed_vf !
!                   ******************************
!
     &(mesh,liebor,ksort,elay,v2dpar,qsx,qsy,ava,npoin,nseg,nptfr,
     & dt,qs,t1,t2,t3,breach,csf_sable,nubo,vnoin)
!
!***********************************************************************
! SISYPHE   V7P0                                   03/06/2014
!***********************************************************************
!
!brief    NON ERODABLE METHOD FOR FINITE VOLUMES.
!
!history  M. GONZALES DE LINARES
!+        07/05/2002
!+        V5P3
!+  First version.
!
!history  J-M HERVOUET (EDF,LNHE)
!+        31/01/2008
!+        V6P0
!+   CORRECTED INITIALISATION ERROR FOR T1 AND T2
!
!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
!+  Names of variables .
!+
!history  R.ATA (EDF-LNHE)
!+        02/06/2014
!+        V7P0
!+  Corrections of normals and nubo tables
!+  after changes in FV data structure of Telemac-2D.
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AVA            |-->| PERCENT AVAILABLE
!| BREACH         |<->| INDICATOR FOR NON ERODIBLE BED (FINITE VOLUMES SCHEMES)
!| DT             |-->| TIME STEP
!| ELAY           |<->| THICKNESS OF SURFACE LAYER
!| KSORT          |-->| CONVENTION FOR FREE OUTPUT
!| LIEBOR         |<->| PHYSICAL BOUNDARY CONDITIONS FOR BED EVOLUTION
!| MESH           |<->| MESH STRUCTURE
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NSEG           |-->| NUMBER OF SEGMENTS PER CONTROL SECTION
!| NUBO           |-->| GLOBAL NUMBER OF EDGE EXTREMITIES
!| QS             |<->| BEDLOAD TRANSPORT RATE
!| QSX            |-->| SOLID DISCHARGE X
!| QSY            |-->| SOLID DISCHARGE Y
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T3             |<->| WORK BIEF_OBJ STRUCTURE
!| V2DPAR         |-->| INTEGRAL OF TEST FUNCTIONS, ASSEMBLED IN PARALLEL
!| VNOIN          |-->| OUTWARD UNIT NORMALS
!| CSF_SABLE      |-->| VOLUME CONCENTRATION OF SAND (1-POROSITY)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe, ex_bedload_nerbed_vf => bedload_nerbed_vf
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_mesh),  INTENT(INOUT) :: MESH
      TYPE(bief_obj),   INTENT(IN)    :: LIEBOR
      TYPE(bief_obj),   INTENT(IN)    :: QSX, QSY
      INTEGER,          INTENT(IN)    :: NPOIN, NSEG, NPTFR,KSORT
      DOUBLE PRECISION, INTENT(IN)    :: DT
      TYPE(bief_obj),   INTENT(INOUT) :: QS, T1, T2, T3
      TYPE(bief_obj),   INTENT(INOUT) :: BREACH
      DOUBLE PRECISION, INTENT(IN)    :: ELAY(npoin),V2DPAR(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: AVA(npoin), CSF_SABLE
! RA
      INTEGER, INTENT(IN)             :: NUBO(2,nseg)
      DOUBLE PRECISION, INTENT(IN)    :: VNOIN(3,nseg)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          :: I, K
      INTEGER          :: IEL, IEL1, IEL2, ISEGIN
      DOUBLE PRECISION :: QSP1, QSP2, QSPC
      DOUBLE PRECISION :: XN, YN, TEMP,PROD_SCAL
      DOUBLE PRECISION :: VNOIN1, VNOIN2, RNORM
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
      ! ****************** !
      ! I - INITIALISES    !
      ! ****************** !
!
      ! BREACH INDICATES IF NON ERODABLE BED WILL BE REACHED
      ! DURING TIME STEP FOR THIS POINT
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     GIVE T1 AND T2 THE SAME STRUCTURE AS QS
!
      CALL cpstvc(qs,t1)
      CALL cpstvc(qs,t2)
!
      DO iel = 1, npoin
        breach%I(iel) = 0
        t1%R(iel)=0.d0
        t2%R(iel)=0.d0
      ENDDO
!
      ! ************************************************* !
      ! II - DETERMINES THE OUTGOING FLUX FOR EACH CELL   ! (_IMP_)
      ! ************************************************* !
      ! THE PRINCIPLE IS THAT QS IS CALCULATED FOR EACH SEGMENT AS
      ! HALF THE SUM OF THE QS AT THE CENTERS OF THE ELEMENTS WHICH
      ! SEGMENT FORMS THE BOUNDARY. IT IS THEN PROJECTED ON THE NORMAL
      ! TO THE SEGMENT, MULTIPLIED BY THE LENGTH OF THE SEGMENT, AND
      ! THIS FLUX TERM IS ADDED (OR SUBTRACTED) TO THE TWO ELEMENTS.
!
      DO isegin = 1, nseg
!
        iel1 = nubo(1,isegin)
        iel2 = nubo(2,isegin)
!
        ! II.1 - SEGMENT LENGTH (RNORM)
        ! ----------------------------------
        vnoin1 = vnoin(1,isegin)
        vnoin2 = vnoin(2,isegin)
        rnorm  = vnoin(3,isegin)
        prod_scal= (mesh%X%R(iel2)-mesh%X%R(iel1))*vnoin1+
     &             (mesh%Y%R(iel2)-mesh%Y%R(iel1))*vnoin2
        IF(prod_scal.LT.0.d0)THEN
          iel1 = nubo(2,isegin)
          iel2 = nubo(1,isegin)
        ENDIF
!
        ! II.2 - PROJECTS QS FOR THE SEGMENT ONTO THE SEGMENT NORMAL
        ! ------------------------------------------------------------
        qsp1 = vnoin1*qsx%R(iel1) + vnoin2*qsy%R(iel1)
        qsp2 = vnoin1*qsx%R(iel2) + vnoin2*qsy%R(iel2)
        qspc = (qsp1+qsp2)*0.5d0
!
        ! II.3 - QS SUCH AS THE OUTGOING FLUX IS MAXIMUM
        ! ----------------------------------------------
        t1%R(iel1) = t1%R(iel1) + rnorm*max(qspc,qsp1,0.d0)
        t1%R(iel2) = t1%R(iel2) - rnorm*min(qspc,qsp2,0.d0)
!
        IF(qspc > 0.d0) THEN
          t2%R(iel1) = t2%R(iel1) + rnorm*qsp1
        ELSEIF(qspc < 0.d0) THEN
          t2%R(iel2) = t2%R(iel2) - rnorm*qsp2
        ENDIF
!
      ENDDO
!
      ! ************************************** !
      ! III - LOOP ON THE BOUNDARY NODES       !
      ! ************************************** !
!
      DO k = 1, nptfr
        iel = mesh%NBOR%I(k)
!
        ! III.1 - FREE EVOLUTION: SEDIMENTS ARE FREE TO LEAVE
        ! ---------------------------------------------------------
        IF (liebor%I(k) == ksort) THEN
!
          ! XNEBOR (*+NPTFR) AND YNEBOR (*+NPTFR)
          ! CONTAIN THE VECTOR NORMAL TO A BOUNDARY NODE
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          xn   = mesh%XNEBOR%R(k+nptfr)
          yn   = mesh%YNEBOR%R(k+nptfr)
          temp = qsx%R(iel)*xn + qsy%R(iel)*yn
          IF (temp > 0.d0) THEN
            t1%R(iel) = t1%R(iel) + temp
            t2%R(iel) = t2%R(iel) + temp
          ENDIF
!
        ENDIF
!
        ! III.2 - FOR A SOLID BOUNDARY: NOTHING TO PROGRAM
        !         BECAUSE THE SEDIMENT FLUX IS ZERO HERE
        !         FOR IMPOSED EVOLUTION : SEE BEDLOAD_SOLVS_VF.F
        ! --------------------------------------------------------
      ENDDO
!
      IF(ncsize > 1) THEN
        CALL parcom(t1, 2, mesh)
        CALL parcom(t2, 2, mesh)
      ENDIF
!
      ! ************************************************ !
      ! IV - COMPUTES THE MAXIMUM FLUX AUTHORISED PER CELL!
      ! ************************************************ !
!
      DO i = 1, npoin
!
        t3%R(i)=elay(i)*v2dpar(i)*ava(i)* csf_sable/dt
        IF (t3%R(i) < 0.d0) t3%R(i) = 0.d0
!
        ! IF THE OUTGOING FLUX IS TOO LARGE, QS IS CAPPED AT THE NODE
        ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        IF(t1%R(i) > t3%R(i)) THEN
          breach%I(i) = 1
          IF(t2%R(i) > t3%R(i)) THEN
            qs%R(i) = qs%R(i)*t3%R(i)/t2%R(i)
          ENDIF
        ENDIF
!
      ENDDO
!
!======================================================================!
!======================================================================!
!
      RETURN
      END