positive_depths_nerd.f

Go to the documentation of this file.

!                   *******************************
                    SUBROUTINE positive_depths_nerd
!                   *******************************
!
     &(t1,t2,t4,h,hn,mesh,flodel,compute_flodel,flbor,dt,
     & unsv2d,npoin,gloseg1,gloseg2,nbor,nptfr,
     & smh,yasmh,pluie,rain,optsou,flulim,limpro,hbor,kdir,info,
     & flopoint,namecode,nitmax,makeflulim)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    Suppresses negative depths by a limitation of fluxes, with the
!+        NERD technique.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        19/11/2016
!+        V7P3
!+   First version, taken out of previous positive_depths.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| COMPUTE_FLODEL |-->| IF YES, COMPUTE FLODEL HERE
!| DOFLULIM       |-->| OPTIONAL, IF YES DOES ARRAY FLULIM
!| DT             |-->| TIME STEP
!| FLBOR          |<->| BOUNDARY FLUXES
!| FLODEL         |<->| FLUXES GIVEN BY SEGMENT
!|                |   | MAY BE COMPUTED HERE (SEE COMPUTE-FLODEL)
!|                |   | OR SIMPLY GIVEN. AT THE EXIT, THE REAL FLUX
!|                |   | TRANSMITTED IS GIVEN BACK.
!| FLOPOINT       |-->| FLUXES GIVEN BY POINTS (ELEMENT BY ELEMENT)
!| FLULIM         |<--| PER SEGMENT: PERCENTAGE OF FLUX THAT HAS NOT
!|                |   | BEEN TRANSMITTED AT THE END OF THE ALGORITHM
!| GLOSEG1        |-->| FIRST POINT OF SEGMENTS
!| GLOSEG2        |-->| SECOND POINT OF SEGMENTS
!| H              |<->| NEW DEPTH
!| HBOR           |-->| PRESCRIBED DEPTHS AT BOUNDARIES
!| HN             |-->| OLD DEPTH
!| INFO           |-->| IF YES, PRINTING INFORMATION ON LISTING
!| KDIR           |-->| CONVENTION FOR DIRICHLET BOUNDARY CONDITION
!| LIMPRO         |-->| TYPE OF BOUNDARY CONDITIONS
!|                |   | IF EQUAL TO KDIR: PRESCRIBED DEPTH.
!| MESH           |<->| MESH STRUCTURE
!| NAMECODE       |-->| NAME OF CALLING CODE (SISYPHE, TELEMEC2D, ETC.)
!| NBOR           |-->| GLOBAL NUMBERS OF BOUNDARY POINTS
!| NITMAX         |-->| MAXIMUM NUMBER OF ITERATIONS
!| NPOIN          |-->| NUMBER OF POINTS IN THE MESH
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| OPTSOU         |-->| OPTION FOR SOURCES 1: NORMAL 2: DIRAC
!| PLUIE          |-->| RAIN IN A BIEF_OBJ, IN M/S.
!| RAIN           |-->| IF YES, THERE IS RAIN OR EVAPORATION
!| SMH            |-->| SOURCE TERMS
!| T1             |-->| WORK ARRAY
!| T2             |-->| WORK ARRAY
!| T4             |-->| WORK ARRAY
!| UNSV2D         |-->| INVERSE OF INTEGRAL OF BASIS FUNCTIONS
!| YASMH          |-->| IF(YES) SMH MUST BE TAKEN INTO ACCOUNT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_positive_depths_nerd => positive_depths_nerd
      USE declarations_telemac, ONLY : deja_pdept_nerd,indic_pdept_nerd
      USE interface_parallel
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NPOIN,NPTFR,OPTSOU,KDIR
      INTEGER, INTENT(IN)             :: NITMAX
      INTEGER, INTENT(IN)             :: GLOSEG1(*),GLOSEG2(*)
      INTEGER, INTENT(IN)             :: NBOR(nptfr)
      INTEGER, INTENT(IN)             :: LIMPRO(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DT,HBOR(nptfr)
      DOUBLE PRECISION, INTENT(INOUT) :: FLULIM(*)
      TYPE(bief_mesh),INTENT(INOUT)   :: MESH
      DOUBLE PRECISION, INTENT(INOUT) :: FLOPOINT(mesh%nelem,3)
      TYPE(bief_obj), INTENT(INOUT)   :: T1,T2,T4,FLODEL,H,FLBOR
      TYPE(bief_obj), INTENT(INOUT)   :: PLUIE
      TYPE(bief_obj), INTENT(IN)      :: UNSV2D,HN,SMH
      LOGICAL, INTENT(IN)             :: YASMH,INFO,RAIN
      LOGICAL, INTENT(IN)             :: COMPUTE_FLODEL
      CHARACTER(LEN=24)               :: NAMECODE
      LOGICAL, INTENT(IN)             :: MAKEFLULIM
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,I1,I2,IPTFR,REMAIN,NEWREMAIN,IR,NITER
      INTEGER NELEM,NELMAX,NSEG
      DOUBLE PRECISION C,CPREV,CINIT,HFL1
      DOUBLE PRECISION SURDT,HSEG1,HSEG2,TET,HFL2
      DOUBLE PRECISION, PARAMETER :: TIERS=1.d0/3.d0
!
      DOUBLE PRECISION, PARAMETER :: EPS_FLUX = 1.d-15
      LOGICAL, PARAMETER :: TESTING = .false.
!
!-----------------------------------------------------------------------
!
!     INDIC_PDEPT_NERD WILL BE A LIST OF SEGMENTS WITH NON ZERO FLUXES
!     HSEG IS THE DEPTH SHARED BETWEEN SEGMENTS
!
!-----------------------------------------------------------------------
!
      surdt=1.d0/dt
      nelem=mesh%NELEM
      nelmax=mesh%NELMAX
      nseg=mesh%NSEG
!
      IF(.NOT.deja_pdept_nerd) THEN
        ALLOCATE(indic_pdept_nerd(nseg))
        deja_pdept_nerd=.true.
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(testing) THEN
        c=1.d99
        cinit=1.d99
        DO i=1,npoin
          c    =min(c    ,h%R(i))
          cinit=min(cinit,hn%R(i))
        ENDDO
        IF(ncsize.GT.1) THEN
          c=p_min(c)
          cinit=p_min(cinit)
        ENDIF
        WRITE(lu,*) 'AVANT TRAITEMENT HAUTEURS NEGATIVES, H MIN=',c
        WRITE(lu,*) 'AVANT TRAITEMENT HAUTEURS NEGATIVES, HN MIN=',cinit
        c=0.d0
        cinit=0.d0
        IF(ncsize.GT.1) THEN
          DO i=1,npoin
            c    =c    +h%R(i) *mesh%IFAC%I(i)/unsv2d%R(i)
            cinit=cinit+hn%R(i)*mesh%IFAC%I(i)/unsv2d%R(i)
          ENDDO
          c    =p_sum(c    )
          cinit=p_sum(cinit)
        ELSE
          DO i=1,npoin
            c    =c    +h%R(i) /unsv2d%R(i)
            cinit=cinit+hn%R(i)/unsv2d%R(i)
          ENDDO
        ENDIF
        WRITE(lu,*) 'AVANT TRAITEMENT MASSE INITIALE=',cinit
        WRITE(lu,*) 'AVANT TRAITEMENT MASSE FINALE  =',c
      ENDIF
!
!     CALCUL DES FLUX PAR SEGMENT (T1 SUIVI DE FALSE NON UTILISE)
!     WHEN RECEIVED FLODEL IS NOT ASSEMBLED IN //
!
      IF(compute_flodel) THEN
!       SO FAR HARDCODED OPTION 2
        CALL flux_ef_vf(flodel%R,flopoint,nseg,nelem,nelmax,
     &                  mesh%ELTSEG%I,mesh%ORISEG%I,
     &                  mesh%IKLE%I,.true.,2)
      ENDIF
!
!     ASSEMBLING FLODEL IN PARALLEL
!
      IF(ncsize.GT.1) THEN
        CALL parcom2_seg(flodel%R,flodel%R,flodel%R,
     &                   nseg,1,2,1,mesh,1,11)
      ENDIF
!
!     ON INTERFACE SEGMENTS, ONLY ONE OF THE TWO TWIN SEGMENTS
!     WILL RECEIVE THE TOTAL FLUX, THE OTHER WILL GET 0.
!
      IF(ncsize.GT.1) THEN
        CALL mult_interface_seg(flodel%R,mesh%NH_COM_SEG%I,
     &                          mesh%NH_COM_SEG%DIM1,
     &                          mesh%NB_NEIGHB_SEG,
     &                          mesh%NB_NEIGHB_PT_SEG%I,
     &                          mesh%LIST_SEND_SEG%I,nseg)
      ENDIF
!
      CALL cpstvc(h,t2)
      CALL cpstvc(h,t4)
      CALL cpstvc(h,t1)
!
      IF(makeflulim) THEN
        DO i=1,nseg
!         SAVING INITIAL FLODEL INTO FLULIM
          flulim(i)=flodel%R(i)
        ENDDO
      ENDIF
!
!     ADDING THE SOURCES (SMH IS NATURALLY ASSEMBLED IN //)
!     FIRST THE POSITIVE SOURCES
!
      IF(yasmh) THEN
        IF(optsou.EQ.1) THEN
          DO i=1,npoin
            h%R(i)=hn%R(i)+dt*max(smh%R(i),0.d0)
          ENDDO
        ELSEIF(optsou.EQ.2) THEN
          DO i=1,npoin
            h%R(i)=hn%R(i)+dt*max(smh%R(i),0.d0)*unsv2d%R(i)
          ENDDO
        ENDIF
      ELSE
        DO i=1,npoin
          h%R(i)=hn%R(i)
        ENDDO
      ENDIF
!
!     RAIN (POSITIVE PLUIE)
!
      IF(rain) THEN
        DO i=1,npoin
          h%R(i)=h%R(i)+dt*max(pluie%R(i),0.d0)
        ENDDO
      ENDIF
!
!     BOUNDARY FLUXES : ADDING THE ENTERING (NEGATIVE) FLUXES
!     FIRST PUTTING FLBOR (BOUNDARY) IN T2 (DOMAIN)
!
!     T2 WILL BE THE ASSEMBLED FLBOR, INITIALISATION HERE
!     IS USELESS EXCEPT THAT PARCOM MAY ADD UNDEFINED
!     NUMBERS (THAT WILL NOT BE USED BUT THAT WILL STOP
!     A COMPILER... TOO BAD!)
      IF(ncsize.GT.1) CALL os('X=0     ',x=t2)
      CALL osdb( 'X=Y     ' ,t2,flbor,flbor,0.d0,mesh)
!     ASSEMBLING T2 (FLBOR IS NOT ASSEMBLED)
      IF(ncsize.GT.1) CALL parcom(t2,2,mesh)
      DO iptfr=1,nptfr
        i=nbor(iptfr)
        h%R(i)=h%R(i)-dt*unsv2d%R(i)*min(t2%R(i),0.d0)
      ENDDO
!
!     FOR OPTIMIZING THE LOOP ON ELEMENTS, ONLY ELEMENTS
!     WITH NON ZERO FLUXES WILL BE CONSIDERED, THIS LIST
!     WILL BE UPDATED. TO START WITH, ALL ELEMENTS IN THE LIST
!
      remain=nseg
!
      DO i=1,remain
        indic_pdept_nerd(i)=i
      ENDDO
!
      cprev=0.d0

!     INITIAL SUM OF FLUXES
      DO i=1,nseg
        cprev=cprev+abs(flodel%R(i))
      ENDDO
      IF(ncsize.GT.1) cprev=p_sum(cprev)
      IF(testing) WRITE(lu,*) 'INITIAL SUM OF FLUXES=',cprev
      cinit=cprev
!
!     LOOP OVER THE LOOP OVER THE ELEMENTS
!
      niter = 0
777   CONTINUE
      niter = niter + 1
!
        IF(niter.EQ.1) THEN
          DO i=1,npoin
            t1%R(i)=0.d0
            t4%R(i)=h%R(i)
          ENDDO
        ELSE
!         NOT ALL THE POINTS NEED TO BE INITIALISED NOW
          DO ir=1,remain
            i=indic_pdept_nerd(ir)
            i1=gloseg1(i)
            i2=gloseg2(i)
            t1%R(i1)=0.d0
            t1%R(i2)=0.d0
!           SAVING THE DEPTH
            t4%R(i1)=h%R(i1)
            t4%R(i2)=h%R(i2)
          ENDDO
!         CANCELLING INTERFACE POINTS (SOME MAY BE ISOLATED IN A SUBDOMAIN
!         AT THE TIP OF AN ACTIVE SEGMENT WHICH IS ELSEWHERE)
          IF(ncsize.GT.1) THEN
            DO iptfr=1,nptir
              i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
              t1%R(i)=0.d0
!             SAVING THE DEPTH
              t4%R(i)=h%R(i)
            ENDDO
          ENDIF
        ENDIF
!
!       COMPUTING THE DEMAND FOR EVERY POINT
!       CANCELLING DEPTHS THAT WILL BE DISTRIBUTED TO ACTIVE SEGMENTS
!       I.E. AS SOON AS THERE IS A DEMAND
!       ANYWAY THEY ARE STORED IN T4 THAT WILL BE USED INSTEAD
!
        DO ir=1,remain
          i=indic_pdept_nerd(ir)
          i1=gloseg1(i)
          i2=gloseg2(i)
          IF(flodel%R(i).GT.eps_flux) THEN
            t1%R(i1)=t1%R(i1)+flodel%R(i)
            h%R(i1)=0.d0
          ELSEIF(flodel%R(i).LT.-eps_flux) THEN
            t1%R(i2)=t1%R(i2)-flodel%R(i)
            h%R(i2)=0.d0
          ENDIF
        ENDDO
!
        IF(ncsize.GT.1) THEN
!         DEMAND ASSEMBLED IN PARALLEL
          CALL parcom(t1,2,mesh)
!         FOR ISOLATED POINTS CONNECTED TO AN ACTIVE
!         SEGMENT THAT IS IN ANOTHER SUBDOMAIN
!         H MUST BE CANCELLED, IF NOT, IT IS SHARED
!         SO THAT IT IS FOUND AGAIN ONCE ASSEMBLED
          DO iptfr=1,nptir
            i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
!           AT THIS LEVEL H IS THE SAME AT INTERFACE POINTS
!           NOW H IS SHARED BETWEEN PROCESSORS TO ANTICIPATE
!           THE FINAL PARALLEL ASSEMBLY
            IF(t1%R(i).GT.eps_flux) THEN
!             POINT THAT WILL GIVE
              h%R(i)=0.d0
            ELSE
!             POINT THAT WILL ONLY RECEIVE
!             IN THIS CASE THEIR DEPTH WILL NOT BE DISTRIBUTED
!             IN THE LOOP ON SEGMENTS, IT IS LEFT UNCHANGED
!             H IS SHARED TO ANTICIPATE THE FURTHER PARCOM
              h%R(i)=h%R(i)*mesh%IFAC%I(i)
            ENDIF
          ENDDO
        ENDIF
!
        c=0.d0
        newremain=0
!
!       TRANSFER OF FLUXES
!
        DO ir=1,remain
          i=indic_pdept_nerd(ir)
          i1=gloseg1(i)
          i2=gloseg2(i)
          IF(flodel%R(i).GT.eps_flux) THEN
!           SHARING ON DEMAND
            hseg1=t4%R(i1)*flodel%R(i)/t1%R(i1)
!           END OF SHARING ON DEMAND
            hfl1= dt*unsv2d%R(i1)*flodel%R(i)
            IF(hfl1.GT.hseg1) THEN
              tet=hseg1/hfl1
              h%R(i2)=h%R(i2)+dt*unsv2d%R(i2)*flodel%R(i)*tet
              flodel%R(i)=flodel%R(i)*(1.d0-tet)
              c=c+flodel%R(i)
              newremain=newremain+1
              indic_pdept_nerd(newremain)=i
            ELSE
              h%R(i1)=h%R(i1)+hseg1-hfl1
              h%R(i2)=h%R(i2)+dt*unsv2d%R(i2)*flodel%R(i)
              flodel%R(i)=0.d0
            ENDIF
          ELSEIF(flodel%R(i).LT.-eps_flux) THEN
!           SHARING ON DEMAND
            hseg2=-t4%R(i2)*flodel%R(i)/t1%R(i2)
!           END OF SHARING ON DEMAND
            hfl2=-dt*unsv2d%R(i2)*flodel%R(i)
            IF(hfl2.GT.hseg2) THEN
              tet=hseg2/hfl2
!             GATHERING DEPTHS
              h%R(i1)=h%R(i1)-dt*unsv2d%R(i1)*flodel%R(i)*tet
              flodel%R(i)=flodel%R(i)*(1.d0-tet)
              c=c-flodel%R(i)
              newremain=newremain+1
              indic_pdept_nerd(newremain)=i
            ELSE
!             GATHERING DEPTHS
              h%R(i1)=h%R(i1)-dt*unsv2d%R(i1)*flodel%R(i)
              h%R(i2)=h%R(i2)+hseg2-hfl2
              flodel%R(i)=0.d0
            ENDIF
          ENDIF
!
        ENDDO
!
      remain=newremain
!     SUMMING THE NEW POSITIVE PARTIAL DEPTHS OF INTERFACE POINTS
      IF(ncsize.GT.1) CALL parcom(h,2,mesh)
      IF(ncsize.GT.1) c=p_sum(c)
      IF(testing) WRITE(lu,*) 'FLUX NON PRIS EN COMPTE=',c
!
!     STOP CRITERION
!
      IF(c.NE.cprev.AND.abs(c-cprev).GT.cinit*1.d-9
     &             .AND.c.NE.0.d0) THEN
        cprev=c
        IF(niter.LT.nitmax) GO TO 777
      ENDIF
!
!     ADDING THE SOURCES (SMH IS NATURALLY ASSEMBLED IN //)
!     NOW THE NEGATIVE SOURCES
!
      IF(yasmh) THEN
        IF(optsou.EQ.1) THEN
          DO i=1,npoin
            h%R(i)=h%R(i)+dt*min(smh%R(i),0.d0)
          ENDDO
        ELSEIF(optsou.EQ.2) THEN
          DO i=1,npoin
            h%R(i)=h%R(i)+dt*min(smh%R(i),0.d0)*unsv2d%R(i)
          ENDDO
        ENDIF
      ENDIF
!
!     EVAPORATION (NEGATIVE PLUIE), WITH A POSSIBLE LIMITATION
!
      IF(rain) THEN
        DO i=1,npoin
          IF(-dt*pluie%R(i).LT.h%R(i)) THEN
            h%R(i)=h%R(i)+dt*min(pluie%R(i),0.d0)
          ELSE
            pluie%R(i)=-h%R(i)/dt
            h%R(i)=0.d0
          ENDIF
        ENDDO
      ENDIF
!
!     BOUNDARY FLUXES : ADDING THE EXITING (POSITIVE) FLUXES
!                       WITH A POSSIBLE LIMITATION
!     IMPORTANT: MUST BE DONE AFTER ALL OTHER SOURCES
!
      DO iptfr=1,nptfr
        i=nbor(iptfr)
!                               T2 = // ASSEMBLED FLBOR
        hfl1=dt*unsv2d%R(i)*max(t2%R(i),0.d0)
        IF(hfl1.GT.h%R(i)) THEN
!         FLBOR ACTUALLY TAKEN INTO ACCOUNT
          flbor%R(iptfr)=flbor%R(iptfr)*h%R(i)/hfl1
          h%R(i)=0.d0
        ELSE
          h%R(i)=h%R(i)-hfl1
        ENDIF
        IF(limpro(iptfr).EQ.kdir) THEN
          IF(hbor(iptfr).LT.0.d0) THEN
            WRITE(lu,*) 'NEGATIVE DEPTH PRESCRIBED ON BOUNDARY'
            WRITE(lu,*) 'CHECK YOUR SPECIFIC SUBROUTINE:'
            IF(namecode(1:7).EQ.'SISYPHE') THEN
              WRITE(lu,*) 'BEDLOAD_SOLVS_FE'
            ELSEIF(namecode(1:4).EQ.'GAIA') THEN
              WRITE(lu,*) 'BEDLOAD_SOLVS_FE'
            ELSEIF(namecode(1:9).EQ.'TELEMAC2D') THEN
              WRITE(lu,*) 'BORD'
            ELSEIF(namecode(1:9).EQ.'TELEMAC3D') THEN
              WRITE(lu,*) 'BORD3D'
            ENDIF
            CALL plante(1)
            stop
          ENDIF
!         HERE V2DPAR WOULD BE MORE CONVENIENT THAN UNSV2D...
          flbor%R(iptfr)=flbor%R(iptfr)
     &                  +(h%R(i)-hbor(iptfr))/(dt*unsv2d%R(i))
          h%R(i)= hbor(iptfr)
        ENDIF
      ENDDO
!
      IF(testing) THEN
        c=1.d99
        DO i=1,npoin
          c=min(c,h%R(i))
        ENDDO
        IF(ncsize.GT.1) c=p_min(c)
        WRITE(lu,*) 'APRES TRAITEMENT HAUTEURS NEGATIVES, HMIN=',c
        c=0.d0
        IF(ncsize.GT.1) THEN
          DO i=1,npoin
            c=c+h%R(i)*mesh%IFAC%I(i)/unsv2d%R(i)
          ENDDO
          c=p_sum(c)
        ELSE
          DO i=1,npoin
            c=c+h%R(i)/unsv2d%R(i)
          ENDDO
        ENDIF
        WRITE(lu,*) 'APRES TRAITEMENT MASSE FINALE =',c
      ENDIF
!
!-----------------------------------------------------------------------
!
!     NOW WE WANT:
!     FLULIM TO BE THE PERCENTAGE OF FLUX THAT HAS BEEN TRANSMITTED
!     FLODEL TO BE THE ACTUAL FLUX THAT HAS BEEN TRANSMITTED
!
!     WE START WITH FLULIM=THE ORIGINAL FLODEL
!
!     WORKING ON ASSEMBLED VALUES TO FIND AN AVERAGE FLULIM
!     THAT WILL BE THE SAME AT INTERFACES
!
!     ASSEMBLING THE REMAINING FLUXES IN FLODEL
!
!
!     NOW WE WANT:
!     FLULIM TO BE THE PERCENTAGE OF FLUX THAT HAS BEEN TRANSMITTED
!     FLODEL TO BE THE ACTUAL FLUX THAT HAS BEEN TRANSMITTED
!
!     WE START WITH FLULIM=THE ORIGINAL FLODEL
!
!     WORKING ON ASSEMBLED VALUES TO FIND AN AVERAGE FLULIM
!     THAT WILL BE THE SAME AT INTERFACES
!
      IF(ncsize.GT.1) THEN
        CALL parcom2_seg(flodel%R,flodel%R,flodel%R,
     &                   nseg,1,2,1,mesh,1,11)
        CALL parcom2_seg(flulim,flulim,flulim,
     &                   nseg,1,2,1,mesh,1,11)
      ENDIF
!
      DO i=1,nseg
!       ACTUAL FLUX TRANSMITTED (=ORIGINAL-REMAINING)
        flodel%R(i)=flulim(i)-flodel%R(i)
!       PERCENTAGE OF ACTUAL FLUX WITH RESPECT TO ORIGINAL FLUX
        IF(abs(flulim(i)).GT.eps_flux) THEN
          flulim(i)=flodel%R(i)/flulim(i)
        ELSE
          flulim(i)=0.d0
        ENDIF
      ENDDO
!
      IF(ncsize.GT.1) THEN
!       SHARING AGAIN FLODEL FOR FURTHER USES
!       ON INTERFACE SEGMENTS, ONLY ONE OF THE TWO TWIN SEGMENTS
!       WILL RECEIVE THE TOTAL FLUX, THE OTHER WILL GET 0.
        CALL mult_interface_seg(flodel%R,mesh%NH_COM_SEG%I,
     &                          mesh%NH_COM_SEG%DIM1,
     &                          mesh%NB_NEIGHB_SEG,
     &                          mesh%NB_NEIGHB_PT_SEG%I,
     &                          mesh%LIST_SEND_SEG%I,nseg)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     CHECKING: COMPARING H(N+1) WITH H RECONSTRUCTED WITH THE FLUXES
!               SOURCES LACKING...
!
      IF(testing) THEN
        DO i=1,npoin
          t1%R(i)=0.d0
        ENDDO
        DO i=1,nseg
          i1=gloseg1(i)
          i2=gloseg2(i)
          t1%R(i1)=t1%R(i1)-dt*unsv2d%R(i1)*flodel%R(i)
          t1%R(i2)=t1%R(i2)+dt*unsv2d%R(i2)*flodel%R(i)
        ENDDO
        DO iptfr=1,nptfr
          i=nbor(iptfr)
          t1%R(i)=t1%R(i)-dt*unsv2d%R(i)*flbor%R(iptfr)
        ENDDO
        IF(ncsize.GT.1) CALL parcom(t1,2,mesh)
        DO i=1,npoin
          t1%R(i)=t1%R(i)+hn%R(i)-h%R(i)
        ENDDO
        WRITE(lu,*) 'ERREUR POSITIVE_DEPTHS_NERD=',p_dots(t1,t1,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(info) THEN
        IF(namecode(1:7).EQ.'SISYPHE'.OR.namecode(1:4).EQ.'GAIA') THEN
          WRITE(lu,102) niter
        ELSE
          WRITE(lu,202) niter
        ENDIF
      ENDIF
!
      IF(niter.EQ.nitmax) THEN
        IF(namecode(1:7).EQ.'SISYPHE'.OR.namecode(1:4).EQ.'GAIA') THEN
          WRITE(lu,103) niter
        ELSE
          WRITE(lu,203) niter
        ENDIF
      ENDIF
!
102   FORMAT(' BEDLOAD EQUATION SOLVED IN ',1i5,' ITERATIONS')
202   FORMAT(' POSITIVE DEPTHS OBTAINED IN ',1i5,' ITERATIONS')
103   FORMAT(' BEDLOAD EQUATION SOLVED IN ',1i5,' ITERATIONS = MAXIMUM',
     &  /,'INCREASE MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES')
203   FORMAT(' POSITIVE DEPTHS SOLVED IN ',1i5,' ITERATIONS = MAXIMUM',
     &  /,'INCREASE MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES')
!
!-----------------------------------------------------------------------
!
      RETURN
      END