bilan_sisyphe.f

Go to the documentation of this file.

!                   ************************
                    SUBROUTINE bilan_sisyphe
!                   ************************
!
     &(e,esomt,t1,vcumu,dt,nptfr,
     & info,zfcl_c,zfcl_s,zfcl_ms,
     & nsicla,voltot,
     & numliq,nfrliq,flbcla,lt,nit,npoin,volu2d,csf_sable,masdep,
     & masdept,charr,susp,slide)
!
!***********************************************************************
! SISYPHE   V6P2                                   21/07/2011
!***********************************************************************
!
!brief    COMPUTES THE MASS BALANCE.
!
!note     T2 IS NOT USED
!
!history  CMGDL
!+
!+        V5P9
!+   CHANGED FOR GRADED SEDIMENT
!
!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)
!+        14/02/2012
!+        V6P2
!+  NSICLM and MAXFRO used instead of 10 and 300. New and compatible
!+  computation: flux given as argument, mass computed differently,
!+  and coefficient CSF_SABLE.
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!history  R.KOPMANN (BAW)
!+        15/02/2019
!+        V7P2
!+   Adding mass balance per class
!+   incorporating nestor volumes in mass balance
!+   Adding mass balance with total volume (initial volume - final volume)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CSF_SABLE      |-->| 1-POROSITY
!| DT             |-->| TIME STEP
!| E              |-->| BED EVOLUTION AT A GIVEN TIME STEP
!| ESOMT          |-->| CUMULATED BED EVOLUTION
!| FLBCLA         |-->| BLOCK OF FLUXES AT BOUNDARY FOR EACH CLASS
!| INFO           |-->| IF YES : INFORMATION IS PRINTED
!| LT             |-->| CURRENT TIME STEP
!| MASDEP         |-->| VOLUME DEPOSITED ON THE BOTTOM FOR EACH CLASS
!|                |   | FROM THE BEGINNING
!| MASDEPT        |-->| VOLUME DEPOSITED ON THE BOTTOM FOR EACH CLASS
!|                |   | FOR THIS TIME STEP
!| MASKEL         |-->| MASKING OF ELEMENTS
!| MESH           |<->| MESH STRUCTURE
!| NFRLIQ         |-->| NUMBER OF LIQUID BOUNDARIES
!| NIT            |-->| NUMBER OF TIME STEPS
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY NODES
!| NSICLA         |-->| NUMBER OF SIZE CLASSES FOR BED MATERIALS
!| NUMLIQ         |-->| LIQUID BOUNDARY NUMBER OF BOUNDARY POINTS
!| QSCLXC         |<->| TRANSPORT RATE FOR EACH CLASS X-DIRECTION
!| QSCLYC         |<->| TRANSPORT RATE FOR EACH CLASS Y-DIRECTION
!| S              |-->| VOID STRUCTURE
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| VCUMU          |<->| VOLUME OF SEDIMENT ENTERING THE DOMAIN
!| VF             |-->| IF YES : FINITE VOLUMES IF NO : FINITE ELEMENTS
!| VOLTOT         |-->| VOLUME TOTAL PER CLASS OF SEDIMENT
!| VOLU2D         |-->| INTEGRAL OF TEST FUNCTIONS (NOT ASSEMBLED IN //)
!| ZFCL_C         |<->| BEDLOAD EVOLUTION FOR EACH SEDIMENT CLASS
!| ZFCL_S         |<->| SUSPENDED LOAD EVOLUTION FOR EACH SEDIMENT CLASS
!| ZFCL_MS        |<->| SLIDE EVOLUTION FOR EACH SEDIMENT CLASS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_sisyphe, ONLY : nsiclm,maxfro
      USE declarations_sisyphe, ONLY : nestor,zf,zr,volini,
     &     vcumucl,rmascl,volnestorcl,volnestorcla
!
      USE declarations_special
      USE interface_parallel, ONLY : p_dsum
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)          :: NPTFR,NFRLIQ,NSICLA,LT,NIT
      INTEGER, INTENT(IN)          :: NPOIN,NUMLIQ(nptfr)
      DOUBLE PRECISION, INTENT(IN) :: DT
      LOGICAL, INTENT(IN)          :: INFO,SUSP,SLIDE,CHARR
!
      DOUBLE PRECISION, INTENT(INOUT) :: VCUMU
      DOUBLE PRECISION, INTENT(IN)    :: CSF_SABLE,VOLTOT(nsicla)
      DOUBLE PRECISION, INTENT(IN)    :: MASDEP(nsicla)
      DOUBLE PRECISION, INTENT(INOUT) :: MASDEPT(nsicla)
!
!-----------------------------------------------------------------------
!
!     VECTOR STRUCTURES
!
      TYPE(bief_obj), INTENT(IN)    :: ZFCL_C
      TYPE(bief_obj), INTENT(IN)    :: E,ESOMT,VOLU2D,ZFCL_S
      TYPE(bief_obj), INTENT(IN)    :: ZFCL_MS
      TYPE(bief_obj), INTENT(INOUT) :: T1,FLBCLA
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,IFRLIQ,IPTFR,ICLA,J
      DOUBLE PRECISION RMASSE,RCUMU,RMASCLA(nsiclm)
      DOUBLE PRECISION VCUMUCLA(nsiclm),FLUXT,FLUXTCLA,VOLDEP
      DOUBLE PRECISION FLT_BOUND(maxfro),VOLDEPC
      DOUBLE PRECISION VOLNESTOR,VOLLOSTPERC,VOLTOT1
      DOUBLE PRECISION VOLINI1,VCUMU1,VOLNESTOR1
!
!-----------------------------------------------------------------------
!
!     COMPUTES THE EVOLUTION (E)
!
      rmasse=0.d0
      DO i=1,npoin
        rmasse=rmasse+e%R(i)*volu2d%R(i)
      ENDDO
      IF(ncsize.GT.1) rmasse = p_dsum(rmasse)
!
!=======================================================================
!
!     COMPUTES THE INTEGRAL OF EVOLUTION AT THE END (ESOMT)
!
      IF(lt.EQ.nit) THEN
        rcumu=0.d0
        DO i=1,npoin
          rcumu=rcumu+esomt%R(i)*volu2d%R(i)
        ENDDO
        IF(ncsize.GT.1) rcumu = p_dsum(rcumu)
      ENDIF
!
!=======================================================================
!
!     COMPUTES THE FLUXES AT THE BOUNDARIES
!
      IF(charr) THEN
        CALL os('X=Y     ',x=t1,y=flbcla%ADR(1)%P)
        IF(nsicla.GT.1) THEN
          DO i=2,nsicla
            CALL os('X=X+Y   ',x=t1,y=flbcla%ADR(i)%P)
          ENDDO
        ENDIF
      ELSE
        CALL cpstvc(flbcla%ADR(1)%P,t1)
        CALL os('X=0     ',x=t1)
      ENDIF
!
      fluxt=0.d0
!
      IF(nfrliq.GT.0) THEN
        DO ifrliq=1,nfrliq
          flt_bound(ifrliq)=0.d0
        ENDDO
        IF(nptfr.GT.0) THEN
          DO iptfr=1,nptfr
            ifrliq=numliq(iptfr)
            IF(ifrliq.GT.0) THEN
              flt_bound(ifrliq)=flt_bound(ifrliq)+t1%R(iptfr)
            ENDIF
          ENDDO
        ENDIF
        IF(ncsize.GT.1) THEN
          DO ifrliq=1,nfrliq
            flt_bound(ifrliq)=p_dsum(flt_bound(ifrliq))
          ENDDO
        ENDIF
        DO ifrliq=1,nfrliq
          fluxt=fluxt+flt_bound(ifrliq)
        ENDDO
      ENDIF
!
      vcumu = vcumu - fluxt*dt/csf_sable
!
!     BALANCE IN EXTENDED GRANULOMETRY
!
      IF(nsicla.GT.1) THEN
!
        DO icla=1,nsicla
!
!       COMPUTES THE EVOLUTION PER CLASS
!
        rmascla(icla)=0.d0
        IF(susp.AND.slide.AND.charr) THEN
          DO i=1,npoin
            rmascla(icla)=rmascla(icla)
     &                   +( zfcl_c%ADR(icla)%P%R(i)
     &                     +zfcl_s%ADR(icla)%P%R(i)
     &                     +zfcl_ms%ADR(icla)%P%R(i) )*volu2d%R(i)
          ENDDO
        ELSEIF(slide.AND.charr) THEN
          DO i=1,npoin
            rmascla(icla)=rmascla(icla)
     &                   +( zfcl_c%ADR(icla)%P%R(i)
     &                     +zfcl_ms%ADR(icla)%P%R(i) )*volu2d%R(i)
          ENDDO
        ELSEIF(susp.AND.charr) THEN
          DO i=1,npoin
            rmascla(icla)=rmascla(icla)
     &                   +( zfcl_c%ADR(icla)%P%R(i)
     &                     +zfcl_s%ADR(icla)%P%R(i) )*volu2d%R(i)
          ENDDO
        ELSEIF(susp.AND.slide) THEN
          DO i=1,npoin
            rmascla(icla)=rmascla(icla)
     &                   +( zfcl_s%ADR(icla)%P%R(i)
     &                     +zfcl_ms%ADR(icla)%P%R(i) )*volu2d%R(i)
          ENDDO
        ELSEIF(susp) THEN
          DO i=1,npoin
            rmascla(icla)=rmascla(icla)
     &                   +zfcl_s%ADR(icla)%P%R(i)*volu2d%R(i)
          ENDDO
        ELSEIF(slide) THEN
          DO i=1,npoin
            rmascla(icla)=rmascla(icla)
     &                   +zfcl_ms%ADR(icla)%P%R(i)*volu2d%R(i)
          ENDDO
        ELSEIF(charr) THEN
          DO i=1,npoin
            rmascla(icla)=rmascla(icla)
     &                   +zfcl_c%ADR(icla)%P%R(i)*volu2d%R(i)
          ENDDO
        ENDIF
        IF(ncsize.GT.1) rmascla(icla) = p_dsum(rmascla(icla))
        rmascl(icla) = rmascl(icla)+rmascla(icla)

!
!       COMPUTES THE FREE FLUXES BY CLASS
!
        fluxtcla=0.d0
        IF(nfrliq.GT.0.AND.charr) THEN
          IF(nptfr.GT.0) THEN
            DO iptfr=1,nptfr
              ifrliq=numliq(iptfr)
              IF(ifrliq.GT.0) THEN
                fluxtcla=fluxtcla+flbcla%ADR(icla)%P%R(iptfr)
              ENDIF
            ENDDO
          ENDIF
          IF(ncsize.GT.1) fluxtcla=p_dsum(fluxtcla)
        ENDIF
!
        vcumucla(icla) = - fluxtcla*dt/csf_sable
        vcumucl(icla) = vcumucl(icla)- fluxtcla*dt/csf_sable
!
        ENDDO
!
      ENDIF
!
!=======================================================================
!
!     GRAIN-FEEDING
!
      voldepc=0.d0
      IF(susp) THEN
        DO i=1,nsicla
          voldepc=voldepc+masdept(i)
        ENDDO
        voldepc=voldepc/csf_sable
      ENDIF
!
!     WRITES OUT THE BALANCE
!
      IF(info) THEN
!
!       GLOBAL BALANCE
!
        WRITE(lu,*)
        WRITE(lu,2000)
        WRITE(lu,2010) rmasse
        IF(nfrliq.GT.0) THEN
          DO ifrliq=1,nfrliq
            WRITE(lu,2110) ifrliq,-flt_bound(ifrliq)/csf_sable
          ENDDO
          WRITE(lu,2111) -fluxt/csf_sable
          IF(susp) WRITE(lu,2112) voldepc
        ENDIF
        volnestor = 0.d0
        IF(nestor) THEN
          DO i=1,nsicla
            volnestor = volnestor + volnestorcla(i)
          END DO
          WRITE(lu,*)'NESTOR VOLUME CHANGE      =   ', volnestor
        ENDIF
        WRITE(lu,2033) rmasse+dt*fluxt/csf_sable-voldepc-volnestor
!       BALANCE PER CLASS
!
        IF(nsicla.GT.1) THEN
          DO i=1,nsicla
            WRITE(lu,*)
            WRITE(lu,*) 'MASS BALANCE FOR SEDIMENT CLASS :',i
            WRITE(lu,*) 'TOTAL VOLUME:',voltot(i)
            WRITE(lu,3010) rmascla(i)
            WRITE(lu,3031) vcumucla(i)
            WRITE(lu,*) 'NESTOR VOLUME PER CLASS: ',volnestorcla(i)
            IF(susp) THEN
              WRITE(lu,3034) masdept(i)/csf_sable
              WRITE(lu,2033) rmascla(i)-vcumucla(i)
     &                                 -masdept(i)/csf_sable
            ELSE
              WRITE(lu,2033) rmascla(i)-vcumucla(i)-volnestorcla(i)
            ENDIF
          ENDDO
        ENDIF
!
!       FINAL GLOBAL BALANCE
!
        IF(lt.EQ.nit) THEN
          WRITE(lu,*)
          WRITE(lu,*)'---------------------------'
          WRITE(lu,*)'FINAL SEDIMENT MASS BALANCE'
          WRITE(lu,*)'---------------------------'
          voldep=0.d0
          volnestor = 0.d0
          DO i=1,nsicla
            voldep=voldep+masdep(i)
          ENDDO
          IF(nestor) THEN
            DO i=1,nsicla
              volnestor=volnestor+volnestorcl(i)
            ENDDO
          ENDIF
          voldep=voldep/csf_sable
          WRITE(lu,*)
          WRITE(lu,2030) rcumu
          WRITE(lu,2031) vcumu
          IF(susp) WRITE(lu,2032) voldep
          IF(nestor) WRITE(lu,*)'CUMULATED NESTOR VOLUME:        :',
     &        volnestor
          WRITE(lu,2033) rcumu-vcumu-voldep-volnestor
          IF(nsicla.EQ.1) THEN
            volini1 = volini(1)
            voltot1 = 0.d0
            DO j=1,npoin
              voltot1 = voltot1 + (zf%R(j)-zr%R(j))*volu2d%R(j)
            ENDDO
            IF(ncsize>1)voltot1 = p_dsum(voltot1)
            volnestor1 = volnestor
            vcumu1 = vcumu
          ELSE ! NSICLA >1
            volini1 = 0.d0
            voltot1 = 0.d0
            volnestor1 = 0.d0
            vcumu1 = 0.d0
            DO i=1,nsicla
              volini1 = volini1+volini(i)
              voltot1 = voltot1+voltot(i)
              volnestor1 = volnestor1+volnestorcl(i)
              vcumu1 = vcumu1+vcumucl(i)
            END DO
          ENDIF

          WRITE(lu,*)'INITIAL VOLUME                  :',volini1
          WRITE(lu,*)'FINAL VOLUME                    :',voltot1
          IF(volini1.GT.0.d0) THEN
            vollostperc = (voltot1-volini1-vcumu1-
     &              volnestor)*100.d0/volini1
          ELSE
            vollostperc = 0.d0
          ENDIF
          WRITE(lu,*)'TOTAL VOLUME LOST               :',
     &        voltot1-volini1-vcumu1-volnestor
          WRITE(lu,*)'                                :'
     &               ,vollostperc,'%'
          IF(nsicla>1) THEN
            DO i=1,nsicla
              WRITE(lu,*)
              WRITE(lu,*) 'MASS BALANCE FOR SEDIMENT CLASS :',i
              WRITE(lu,*)'VOLUME THAT ENTERED THE DOMAIN PER CLASS   ',
     &                    vcumucl(i)
              WRITE(lu,*)'SUM OF THE CUMULATED EVOLUTIONS PER CLASS  ',
     &                    rmascl(i)
              WRITE(lu,*)'BALANCE FLUXES AND EVOLUTIONS PER CLASS    ',
     &                    rmascl(i)-vcumucl(i)
              WRITE(lu,*) ''
              WRITE(lu,*)'INITIAL VOLUME PER CLASS OVER ALL LAYERS   ',
     &                    volini(i)
              WRITE(lu,*)'FINAL VOLUME  PER CLASS OVER ALL LAYERS    ',
     &                    voltot(i)
              IF(nestor)
     &           WRITE(lu,*)
     &          'DREDGED/DISPOSED VOLUME PER CLASS          ',
     &                    volnestorcl(i)
              IF(volini(i).GT.0.d0) THEN
                vollostperc = (voltot(i)-volini(i)-vcumucl(i)-
     &           volnestorcl(i))*100.d0/volini(i)
              ELSE
                vollostperc = 0.d0
              ENDIF
              WRITE(lu,*)'TOTAL VOLUME LOST PER CLASS OVER ALL LAYERS',
     &                    voltot(i)-volini(i)-vcumucl(i)-volnestorcl(i)
              WRITE(lu,*)'                                           '
     &            ,vollostperc,'%'
            END DO
          ENDIF
        ENDIF ! lt==nit
!
!       IF(LGRAFED) THEN
!         WRITE(LU, 4001) MASST
!         WRITE(LU, 4011) MASS_GF
!       ENDIF
!
      ENDIF
!
2000  FORMAT(1x,'MASS-BALANCE (IN VOLUME, INCLUDING VOID): ')
2010  FORMAT(1x,'SUM OF THE EVOLUTIONS : ',g16.7,' M3')
2030  FORMAT(1x,'SUM OF THE CUMULATED EVOLUTIONS : ',g16.7)
2031  FORMAT(1x,'VOLUME THAT ENTERED THE DOMAIN  : ',g16.7,' M3'
     &         ,'  ( IF <0 EXIT )')
2032  FORMAT(1x,'VOLUME DEPOSITED ON THE BOTTOM  : ',g16.7,' M3'
     &         ,'  ( IF <0 ERODED )')
2033  FORMAT(1x,'LOST VOLUME                     : ',g16.7,' M3'
     &         ,'  ( IF <0 EXIT )')
2110  FORMAT(1x,'BOUNDARY ',1i3,' BEDLOAD FLUX = ',g16.7,
     &          '  ( M3/S  >0 = ENTERING )')
2111  FORMAT(1x,'TOTAL        BEDLOAD FLUX = ',g16.7,
     &          '  ( M3/S  >0 = ENTERING )')
2112  FORMAT(1x,'DEPOSIT ON BOTTOM         = ',g16.7,
     &          '  ( M3/S )')
3010  FORMAT(1x,'SUM OF THE EVOLUTIONS FOR THIS CLASS: ',g16.7)
3031  FORMAT(1x,'VOLUME THAT ENTERED THE DOMAIN FOR THIS CLASS: '
     &       ,g16.7,' M3')
3034  FORMAT(1x,'VOLUME DEPOSITED ON BOTTOM FOR THIS CLASS:     '
     &       ,g16.7,' M3')
!
!-----------------------------------------------------------------------
!
      RETURN
      END