bedload_solvs_vf.f

Go to the documentation of this file.

!                   ***************************
                    SUBROUTINE bedload_solvs_vf
!                   ***************************
!
     &(mesh,qsx,qsy,limtec,unsv2d,ebor,breach,nseg,nptfr,npoin,
     & kent,kdir,kddl,dt,zfcl,flux,csf_sable,flbcla,ava,liqbor,qbor,
     & nubo,vnoin)
!
!***********************************************************************
! SISYPHE   V7P0                                      03/06/2014
!***********************************************************************
!
!brief    SOLVES EXNER EQUATION WITH THE FINITE VOLUME METHOD.
!
!history  M. GONZALES DE LINARES
!+        07/05/2002
!+        V5P5
!+   First version.
!
!history  J-M HERVOUET (EDF, LNHE)
!+        30/10/2007
!+        V5P8
!+   UNSV2D +DIRICL DELETED
!
!history  JMH
!+        15/09/2009
!+
!+   KDIR KDDL ADDED (WERE HARD-CODED BEFORE !!!)
!
!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)
!+        21/02/2012
!+        V6P2
!+  Corrections for a perfect mass balance: FLBTRA built to be used in
!+  bilan_sisyphe, coefficient AVA added in Dirichlet value, QBOR
!+  dealt with.
!
!history  R.ATA (EDF-LNHE)
!+        02/06/2014
!+        V7P0
!+  Corrections of normals and nubo tables
!+  after changes in FV data structure of Telemac2d
!+
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AVA            |-->| PERCENTAGE OF CLASS IN SEDIMENT
!| BREACH         |<->| INDICATOR FOR NON ERODIBLE BED
!| DT             |-->| TIME STEP
!| EBOR           |<->| BOUNDARY CONDITION FOR BED EVOLUTION (DIRICHLET)
!| FLBCLA         |<->| FLUXES AT BOUNDARY FOR THE CLASS
!| FLUX           |<->| SEDIMENT FLUX
!| KDIR           |-->| CONVENTION FOR DIRICHLET VALUE
!| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
!| KENT           |-->| CONVENTION FOR PRESCRIBED VALUE AT ENTRANCE
!| LIMTEC         |-->| TECHNICAL BOUNDARY CONDITIONS FOR BED EVOLUTION
!| LIQBOR         |-->| TYPE OF BOUNDARY CONDITIONS FOR BEDLOAD DISCHARGE
!| MESH           |<->| MESH STRUCTURE
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NSEG           |-->| NUMBER OF SEGMENTS
!| NUBO           |-->| GLOBAL NUMBER OF EDGE EXTREMITIES
!| QBOR           |-->| PRESCRIBED BEDLOAD DISCHARGES
!| QSX            |<->| BEDLOAD TRANSPORT RATE X-DIRECTION
!| QSY            |<->| BEDLOAD TRANSPORT RATE Y-DIRECTION
!| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
!| VNOIN          |-->| OUTWARD UNIT NORMALS
!| ZFCL           |<->| BEDLOAD EVOLUTION FOR EACH SEDIMENT CLASS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe, ex_bedload_solvs_vf => bedload_solvs_vf
      USE bief
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_mesh),  INTENT(INOUT) :: MESH
      TYPE(bief_obj),   INTENT(IN)    :: QSX,QSY,LIMTEC,UNSV2D,EBOR
      TYPE(bief_obj),   INTENT(IN)    :: BREACH,LIQBOR,QBOR
      INTEGER,          INTENT(IN)    :: NSEG,NPTFR,NPOIN,KENT,KDIR,KDDL
      DOUBLE PRECISION, INTENT(IN)    :: DT,CSF_SABLE
      DOUBLE PRECISION, INTENT(IN)    :: AVA(npoin)
      TYPE(bief_obj),   INTENT(INOUT) :: FLBCLA,ZFCL,FLUX
      INTEGER, INTENT(IN)             :: NUBO(2,nseg)
      DOUBLE PRECISION, INTENT(IN)    :: VNOIN(3,nseg)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          :: ISEGIN,K,N,IEL,IEL1,IEL2,IELEM,NELEM,ERR
      DOUBLE PRECISION :: QSMOY1,QSMOY2,QSP,VNOIN1,VNOIN2,RNORM,XN,YN
      DOUBLE PRECISION :: ZFCLDIR,PROD_SCAL
      LOGICAL, ALLOCATABLE :: YESNO(:)
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
      nelem = mesh%NELEM
      ALLOCATE(yesno(nseg),stat=err)
!     INITIALIZATION OF YESNO
      DO k=1,nseg
        yesno(k)=.false.
      ENDDO
!
!     INITIALISES THE DIVERGENCE
!
      CALL os('X=0     ', x=flux)
!
!     DETERMINES THE OUTGOING FLUX FOR EACH CELL
!
      DO ielem=1,nelem
        DO k=1,3
          isegin = mesh%ELTSEG%I(ielem+(k-1)*nelem)
          IF(.NOT.yesno(isegin)) THEN
            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 - QS FOR THE SEGMENT, BROKEN UP ACCORDING TO X AND Y
            ! ---------------------------------------------------------
            qsmoy1 = 0.5d0*(qsx%R(iel1) + qsx%R(iel2))
            qsmoy2 = 0.5d0*(qsy%R(iel1) + qsy%R(iel2))
            ! II.3 - PROJECTS QS FOR THE SEGMENT ONTO THE SEGMENT NORMAL
            ! ----------------------------------------------------------
            qsp = vnoin1*qsmoy1 + vnoin2*qsmoy2
            ! II.4 - UPWIND SCHEME ON NODES WITH A "PROBLEM"
            ! ----------------------------------------------
            IF(breach%I(iel1).EQ.1.AND.qsp.GT.0.d0) THEN
              qsmoy1 = qsx%R(iel1)
              qsmoy2 = qsy%R(iel1)
            ENDIF
            IF(breach%I(iel2).EQ.1.AND.qsp.LT.0.d0) THEN
              qsmoy1 = qsx%R(iel2)
              qsmoy2 = qsy%R(iel2)
            ENDIF
            qsp = vnoin1*qsmoy1 + vnoin2*qsmoy2
            ! II.5 - INTEGRATES BY THE SEGMENT LENGTH
            ! ---------------------------------------
            IF(mesh%IFABOR%I(ielem+(k-1)*nelem).NE.-2) THEN
              flux%R(iel1) = flux%R(iel1) + rnorm*qsp
              flux%R(iel2) = flux%R(iel2) - rnorm*qsp
            ELSE
!           IF IT IS AN INTERFACE EDGE IN PARALLEL, ONLY HALF
              flux%R(iel1) = flux%R(iel1) + 0.5d0*rnorm*qsp
              flux%R(iel2) = flux%R(iel2) - 0.5d0*rnorm*qsp
            ENDIF
            yesno(isegin)=.true.
          ENDIF
        ENDDO
      ENDDO
!
!     BOUNDARIES
!
      IF(ncsize.LE.1) THEN
        DO k = 1 , nptfr
          iel = mesh%NBOR%I(k)
!         VECTOR NORMAL TO A BOUNDARY NODE
!         VERSION WHICH IS NOT NORMED
          xn = mesh%XNEBOR%R(k+nptfr)
          yn = mesh%YNEBOR%R(k+nptfr)
!
!         ADDING BOUNDARY FLUX AT OPEN BOUNDARIES
!         QBOR HAS PRIORITY HERE, SO IN CASE OF LIQBOR=KENT
!         LIMTEC=KDIR WILL NOT BE CONSIDERED, SEE BEDLOAD_SOLVS_FE
!         FOR MORE EXPLANATIONS
!
          IF(liqbor%I(k).EQ.kent) THEN
            flbcla%R(k) = qbor%R(k)
            flux%R(iel) = flux%R(iel) + flbcla%R(k)
          ELSEIF(limtec%I(k).EQ.kddl.OR.limtec%I(k).EQ.kdir) THEN
!           IF KDIR WILL BE UPDATED LATER
            flbcla%R(k)= qsx%R(iel)*xn + qsy%R(iel)*yn
!           ADDS THE CONTRIBUTION OF THE FLUX ON THE BOUNDARY SEGMENT
            flux%R(iel) = flux%R(iel) + flbcla%R(k)
          ELSE
!           NO SEDIMENT FLUX ACCROSS SOLID BOUNDARIES
            flbcla%R(k)=0.d0
          ENDIF
        ENDDO
      ELSE
!       COMPUTED ONLY ONCE WHEN THE LIQUID BOUNDARY NODE BELONGS TO
!       2 SUBDDOMAINS
        DO k = 1 , nptfr
          iel = mesh%NBOR%I(k)
!         VECTOR NORMAL TO A BOUNDARY NODE
!         VERSION WHICH IS NOT NORMED
          xn = mesh%XNEBOR%R(k+nptfr)
          yn = mesh%YNEBOR%R(k+nptfr)
!
!         ADDING BOUNDARY FLUX AT OPEN BOUNDARIES
!         QBOR HAS PRIORITY HERE, SO IN CASE OF LIQBOR=KENT
!         LIMTEC=KDIR WILL NOT BE CONSIDERED, SEE BEDLOAD_SOLVS_FE
!         FOR MORE EXPLANATIONS
!
          IF(liqbor%I(k).EQ.kent) THEN
            flbcla%R(k) = qbor%R(k)
            flux%R(iel) = flux%R(iel) + flbcla%R(k)*mesh%IFAC%I(iel)
          ELSEIF(limtec%I(k).EQ.kddl.OR.limtec%I(k).EQ.kdir) THEN
!           IF KDIR WILL BE UPDATED LATER
            flbcla%R(k)= qsx%R(iel)*xn + qsy%R(iel)*yn
!           ADDS THE CONTRIBUTION OF THE FLUX ON THE BOUNDARY SEGMENT
            flux%R(iel) = flux%R(iel) + flbcla%R(k)*mesh%IFAC%I(iel)
          ELSE
!           NO SEDIMENT FLUX ACCROSS SOLID BOUNDARIES
            flbcla%R(k)=0.d0
          ENDIF
        ENDDO
      ENDIF
!
!     ASSEMBLING IN PARALLEL
!
      IF(ncsize.GT.1) CALL parcom(flux, 2, mesh)
!
!     SOLVING, NEGATIVE SIGN BECAUSE OUTGOING FLUX IS POSITIVE
!     TODO: NOTE JMH: FLUX MUST BE HERE DIV(QS)
!
      CALL os('X=CYZ   ', x=zfcl, y=flux, z=unsv2d, c=-dt)
!
      DO k=1,nptfr
!                                  PRIORITY OF LIQBOR, SEE ABOVE
        IF(limtec%I(k).EQ.kdir.AND.liqbor%I(k).NE.kent) THEN
          n = mesh%NBOR%I(k)
!         ZFCLDIR: DIRICHLET VALUE OF EVOLUTION, ZFCL WILL BE DIVIDED BY
!         CSF_SABLE AFTER, AND THEN IT WILL BE AVA(N)*EBOR...
          zfcldir = ava(n)*ebor%R(k)*csf_sable
!         CORRECTION OF BOUNDARY FLUX TO GET ZFCLDIR
          flbcla%R(k)=flbcla%R(k)-(zfcldir-zfcl%R(n))/(dt*unsv2d%R(n))
!         ZFCLDIR FINALLY IMPOSED
          zfcl%R(n) = zfcldir
        ENDIF
!
      ENDDO
!
      DEALLOCATE(yesno)
!
!======================================================================!
!======================================================================!
!
      RETURN
      END