cvtrvf_nerd.f

Go to the documentation of this file.

!                   **********************
                    SUBROUTINE cvtrvf_nerd
!                   **********************
!
     &(f,fn,fscexp,h,hn,hprop,udel,vdel,dm1,zconv,solsys,
     & sm,smh,yasmh,smi,yasmi,fbor,masktr,mesh,
     & t1,t2,t3,t4,t5,t6,t7,ht,dt,entet,
     & msk,maskel,optsou,limtra,kdir,kddl,nptfr,flbor,
     & yaflbor,unsv2d,iopt,flbortra,gloseg1,gloseg2,nbor,
     & flulim,yaflulim,rain,pluie,train,given_flux,flux_given,
     & nitmax)
!
!***********************************************************************
! BIEF   V7P2
!***********************************************************************
!
!brief    NERD advection scheme in 2D.
!
!warning  SEE BELOW FOR DEFINITION OF IOPT1 AND IOPT2, RETRIEVED FROM IOPT
!+        IOPT2=1 NOT TREATED HERE, MASS-CONSERVATION WILL BE DOWNGRADED
!+        IN THIS CASE (A CORRECT TREATMENT MAY RESULT IN INFINITE F)
!+        THE PROGRAM WILL NOT STOP IF IOPT2=1
!
!history  J-M HERVOUET (LNHE)
!+        16/07/2016
!+        V7P2
!+   First version. Actually the former cvtrvf_pos renamed cvtrvf_nerd.
!+   All the part with OPTION = 1 is removed and is now in cvtrvf_eria.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DIFT           |-->| LOGICAL, IF YES THERE IS DIFFUSION OF F
!| DM1            |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
!| DT             |-->| TIME STEP
!| ENTET          |-->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS
!|                |   | CONSERVATION.
!| F              |<--| F AT TIME T(N+1)
!| FBOR           |-->| DIRICHLET CONDITIONS ON F.
!| FLBOR          |-->| FLUXES AT BOUNDARIES
!| FLBORTRA       |<->| TRACER FLUXES AT BOUNDARIES
!| FLUX_GIVEN     |-->| IF GIVEN_FLUX=YES, THE FLUX IS GIVEN IN
!|                |   | GIVEN_FLUX
!| FN             |-->| F AT TIME T(N)
!| FSCEXP         |-->| EXPLICIT PART OF THE SOURCE TERM
!|                |   | EQUAL TO ZERO EVERYWHERE BUT ON SOURCES
!|                |   | WHERE THERE IS FSCE - (1-TETAT) FN
!|                |   | SEE DIFSOU
!| GIVEN_FLUX     |-->| IF GIVEN_FLUX=YES, THE FLUX IS GIVEN IN
!|                |   | GIVEN_FLUX AND WILL NOT BE COMPUTED HERE
!| GLOSEG1        |-->| FIRST POINT OF SEGMENTS
!| GLOSEG2        |-->| SECOND POINT OF SEGMENTS
!| HT             |<--| WORK ARRAY (DEPTH GOING FROM HN TO H)
!| HPROP          |-->| PROPAGATION DEPTH (DONE IN CVDFTR).
!| IOPT           |-->| OPTIONS FOR COMPUTATION (NUMBER BETWEEN 0 AND 13)
!|                |   | THE TENS (IOPT2, I.E. 0 OR 1):
!|                |   | 0: UCONV OBEYS THE CONTINUITY EQUATION
!|                |   | 1: UCONV DOES NOT OBEY THE CONTINUITY EQUATION
!|                |   | THE UNITS (IOPT1, I.E. 0 TO 3): VARIANT FOR FLUXES
!|                |   | 0: CONSTANT PER ELEMENT = 0
!|                |   | 1: CHI-TUAN PHAM'S CONSTANT
!|                |   | 2: N SCHEME
!|                |   | 3: PSI SCHEME
!| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
!| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
!| LIMTRA         |-->| BOUNDARY CONDITIONS ON BOOUNDARY POINTS
!| MASKEL         |-->| MASKING OF ELEMENTS
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
!| NBOR           |-->| GLOBAL NUMBERS OF BOUNDARY POINTS
!| NITMAX         |-->| MAXIMUM NUMBER OF ITERATIONS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| OPTSOU         |-->| TYPE OF SOURCES
!|                |   | 1: NORMAL
!|                |   | 2: DIRAC
!| PLUIE          |-->| RAIN OR EVAPORATION, IN M/S
!| RAIN           |-->| IF YES: RAIN OR EVAPORATION
!| SM             |-->| SOURCE TERMS.
!| SMH            |-->| SOURCE TERM IN CONTINUITY EQUATION
!| SMI            |-->| IMPLICIT SOURCE TERM
!| SOLSYS         |-->| 1 OR 2. IF 2 ADVECTION FIELD IS UCONV + DM1*GRAD(ZCONV)
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T3             |<->| WORK BIEF_OBJ STRUCTURE
!| T4             |<->| WORK BIEF_OBJ STRUCTURE
!| T5             |<->| WORK BIEF_OBJ STRUCTURE
!| T6             |<->| WORK BIEF_OBJ STRUCTURE
!| T7             |<->| WORK BIEF_OBJ STRUCTURE
!| TRAIN          |-->| VALUE OF TRACER IN THE RAIN
!| UDEL           |-->| X-COMPONENT OF ADVECTION VELOCITY
!| UNSV2D         |-->| INVERSE OF V2DPAR
!| VDEL           |-->| X-COMPONENT OF ADVECTION VELOCITY
!| YAFLBOR        |-->| IF YES FLBOR IS GIVEN
!| YASMH          |-->| IF YES, SMH MUST BE TAKEN INTO ACCOUNT
!| YASMI          |-->| IF YES, SMI MUST BE TAKEN INTO ACCOUNT
!| ZCONV          |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_cvtrvf_nerd => cvtrvf_nerd
      USE declarations_telemac, ONLY : deja_cpos, indic_cpos
!
      USE declarations_special
      USE interface_parallel
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: OPTSOU,KDIR,NPTFR,SOLSYS
      INTEGER, INTENT(IN)             :: KDDL,IOPT,NITMAX
      INTEGER, INTENT(IN)             :: GLOSEG1(*),GLOSEG2(*)
      INTEGER, INTENT(IN)             :: NBOR(nptfr)
      INTEGER, INTENT(INOUT)          :: LIMTRA(nptfr)
!                                                                NSEG
      DOUBLE PRECISION, INTENT(IN)    :: DT,TRAIN,FLULIM(*)
      LOGICAL, INTENT(IN)             :: YASMH,YAFLBOR,RAIN
      LOGICAL, INTENT(IN)             :: MSK,ENTET,YASMI,YAFLULIM
      LOGICAL, INTENT(IN)             :: FLUX_GIVEN
      TYPE(bief_obj), INTENT(IN)      :: MASKEL,H,HN,DM1,ZCONV
      TYPE(bief_obj), INTENT(IN)      :: UNSV2D,HPROP
      TYPE(bief_obj), INTENT(INOUT)   :: F,SM,HT
      TYPE(bief_obj), INTENT(IN)      :: UDEL,VDEL,FN,SMI,SMH
      TYPE(bief_obj), INTENT(INOUT)   :: FLBORTRA,FBOR
      TYPE(bief_obj), INTENT(INOUT)   :: T1,T2,T3,T4,T5,T6,T7
      TYPE(bief_obj), INTENT(IN)      :: FSCEXP,MASKTR
      TYPE(bief_obj), INTENT(INOUT)   :: FLBOR
      TYPE(bief_obj), INTENT(IN)      :: PLUIE,GIVEN_FLUX
      TYPE(bief_mesh)                 :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,IOPT1,IOPT2,NPOIN,IPTFR,I1,I2,NITER,NEWREMAIN
      INTEGER IR,N,NELEM,NELMAX,REMAIN,NSEG
!
!-----------------------------------------------------------------------
!
      DOUBLE PRECISION C,CPREV,CINIT,HFL1,HFL2,TET,HSEG1,HSEG2
      CHARACTER(LEN=16) FORMUL
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: FXMAT
      DOUBLE PRECISION, PARAMETER :: TIERS=1.d0/3.d0
      LOGICAL, PARAMETER :: TESTING = .false.
      DOUBLE PRECISION, PARAMETER :: EPS_FLUX = 1.d-15
!
!-----------------------------------------------------------------------
!
      nelem=mesh%NELEM
      nelmax=mesh%NELMAX
      nseg=mesh%NSEG
      npoin=h%DIM1
!
!     INDIC_CPOS WILL BE A LIST OF SEGMENTS WITH NON ZERO FLUXES
!
      IF(.NOT.deja_cpos) THEN
        ALLOCATE(indic_cpos(mesh%NSEG))
        deja_cpos=.true.
      ENDIF
!
!-----------------------------------------------------------------------
!
!     THERE IS PLENTY OF MEMORY AVAILABLE IN A BIEF_MESH STRUCTURE...
!
      fxmat=>mesh%MSEG%X%R(1:nseg)
!
!-----------------------------------------------------------------------
!
!     EXTRACTING OPTIONS, AND CONTROL
!
      iopt2=iopt/10
      iopt1=iopt-10*iopt2
      IF(iopt1.LT.0.OR.iopt1.GT.3) THEN
        WRITE(lu,*) 'CVTRVF_NERD: OPTION IOPT1 UNKNOWN: ',iopt1
        CALL plante(1)
        stop
      ENDIF
      IF(iopt2.NE.0.AND.iopt2.NE.1) THEN
        WRITE(lu,*) 'CVTRVF_NERD: OPTION IOPT2 UNKNOWN: ',iopt2
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     STARTING AGAIN FROM NON CORRECTED DEPTH
!
      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
      ENDIF
!
      CALL cpstvc(h,t1)
      CALL cpstvc(h,t2)
      CALL cpstvc(f,t3)
      CALL cpstvc(h,t4)
      CALL cpstvc(f,t5)
      CALL cpstvc(f,t6)
      CALL cpstvc(h,t7)
!
!     T2 WILL BE THE ASSEMBLED FLBOR, INITIALISATION TO 0 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)
!
!     OPTIONAL COMPUTATION OF FLUXES AT BOUNDARIES
!
      IF(.NOT.yaflbor) THEN
!       MASK=8 FOR LIQUID BOUNDARIES
        CALL vector(flbor,'=','FLUBDF          ',1,1.d0,
     &              hprop,hprop,hprop,
     &              udel , vdel, vdel,mesh,.true.,masktr%ADR(8)%P)
      ENDIF
!
!     BOUNDARY FLUXES : ADDING THE ENTERING (NEGATIVE) FLUXES
!     FIRST PUTTING FLBOR (BOUNDARY) IN T2 (DOMAIN)
      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)
!     POSSIBLE CORRECTION OF LIMTRA AND FBOR (LIMTRA HAS BEEN DONE BY
!     DIFFIN WITH U.N)
      DO i=1,mesh%NPTFR
        n=mesh%NBOR%I(i)
        IF(limtra(i).EQ.kdir.AND.t2%R(n).GT.0.d0) THEN
!         DIRICHLET DECLARED FOR AN EXIT
          limtra(i)=kddl
        ELSEIF(limtra(i).EQ.kddl.AND.t2%R(n).LT.0.d0) THEN
!         DEGREE OF FREEDOM DECLARED FOR AN ENTRANCE
!         IN THIS CASE THE PRESCRIBED VALUE MAY NOT HAVE BEEN GIVEN
!         FN IS TAKEN HERE
          limtra(i)=kdir
          fbor%R(i)=fn%R(n)
        ELSE
!         INITIALISING FLBORTRA FOR OTHER POINTS
          flbortra%R(i)=0.d0
        ENDIF
      ENDDO
!
!     COMPUTING FLUXES
!
      IF(.NOT.flux_given) THEN
        formul='HUGRADP         '
        IF(solsys.EQ.2) formul(8:8)='2'
!       COMPUTING FLUXES LEAVING NODES
        CALL vector(t1,'=',formul,h%ELM,-1.d0,
     &              hprop,dm1,zconv,udel,vdel,vdel,mesh,msk,maskel)
!                   T1 AS HUGRADP IS NOT USED AS AN ASSEMBLED VECTOR
!                   BUT TO GET THE NON ASSEMBLED FORM MESH%W
!       COMPUTING FLUXES BY SEGMENT (TE1 SUIVI DE FALSE NON UTILISE)
!       FXMAT IS NOT ASSEMBLED IN //
        CALL flux_ef_vf(fxmat,mesh%W%R,mesh%NSEG,nelem,nelmax,
     &                  mesh%ELTSEG%I,mesh%ORISEG%I,
     &                  mesh%IKLE%I,.true.,iopt1)
!       LIMITATION OF FLUXES (IF REQUESTED, E.G. USED BY SISYPHE)
        IF(yaflulim) THEN
          DO i=1,mesh%NSEG
            fxmat(i)=fxmat(i)*flulim(i)
          ENDDO
        ENDIF
      ELSE
        DO i=1,mesh%NSEG
          fxmat(i)=given_flux%R(i)
        ENDDO
      ENDIF
!     INTERFACE SEGMENTS: ONLY ONE OF THE TWINS WILL RECEIVE
!     THE ASSEMBLED FLUX
      IF(ncsize.GT.1) THEN
        CALL parcom2_seg(fxmat,fxmat,fxmat,mesh%NSEG,1,2,1,mesh,
     &                   1,11)
        CALL mult_interface_seg(fxmat,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,mesh%NSEG)
      ENDIF
!
!     INITIALIZING F
!
      CALL os('X=Y     ',x=f,y=fn)
      CALL os('X=Y     ',x=ht,y=hn)
!
!     ADDING THE POSITIVE SOURCES (SMH IS NATURALLY ASSEMBLED IN //)
!
      IF(yasmh) THEN
        IF(optsou.EQ.1) THEN
          DO i=1,npoin
            IF(smh%R(i).GT.0.d0) THEN
              ht%R(i)=ht%R(i)+dt*smh%R(i)
              f%R(i)=f%R(i)+dt/max(ht%R(i),1.d-4)*smh%R(i)*fscexp%R(i)
            ENDIF
          ENDDO
        ELSEIF(optsou.EQ.2) THEN
          DO i=1,npoin
            IF(smh%R(i).GT.0.d0) THEN
              ht%R(i)=ht%R(i)+dt*smh%R(i)*unsv2d%R(i)
              f%R(i)=f%R(i)+dt/max(ht%R(i),1.d-4)*
     &                        unsv2d%R(i)*smh%R(i)*fscexp%R(i)
            ENDIF
          ENDDO
        ENDIF
      ENDIF
!
!     RAIN-EVAPORATION: RAIN FIRST, EVAPORATION IN THE END
!
      IF(rain) THEN
        DO i=1,npoin
          c=max(pluie%R(i),0.d0)
          ht%R(i)=ht%R(i)+dt*c
!                                                VALUE IN RAIN
          f%R(i)=f%R(i)+dt/max(ht%R(i),1.d-4)*c*(train-f%R(i))
        ENDDO
      ENDIF
!
      DO iptfr=1,nptfr
        i=nbor(iptfr)
        ht%R(i)=ht%R(i)-dt*unsv2d%R(i)*min(t2%R(i),0.d0)
!       ENTERING FLUXES OF TRACERS (T2<0)
!       THE FINAL DEPTH IS TAKEN
        IF(limtra(iptfr).EQ.kdir) THEN
          f%R(i)=f%R(i)-dt/max(ht%R(i),1.d-4)*
     &       unsv2d%R(i)*t2%R(i)*(fbor%R(iptfr)-f%R(i))
!       ELSEIF(LIMTRA(IPTFR).EQ.KDDL) THEN
!         NOTHING TO DO
        ENDIF
      ENDDO
!
!     FOR OPTIMIZING THE LOOP ON SEGMENTS, ONLY SEGMENTS
!     WITH NON ZERO FLUXES WILL BE CONSIDERED, THIS LIST
!     WILL BE UPDATED. TO START WITH, ALL FLUXES ASSUMED NON ZERO
!
      remain=nseg
!
      DO i=1,remain
        indic_cpos(i)=i
      ENDDO
!
!     MAXIMUM INITIAL FLUX
!
      cprev=0.d0
      DO i=1,mesh%NSEG
        cprev=cprev+abs(fxmat(i))
      ENDDO
      IF(ncsize.GT.1) cprev=p_sum(cprev)
      IF(testing) WRITE(lu,*) 'INITIAL SUM OF FLUXES=',cprev
      cinit=cprev
!
      niter = 0
777   CONTINUE
      niter = niter + 1
!
      IF(niter.EQ.1) THEN
        DO i=1,npoin
          t1%R(i)=0.d0
          t4%R(i)=ht%R(i)
          t6%R(i)=f%R(i)
          t5%R(i)=ht%R(i)*f%R(i)
        ENDDO
        IF(ncsize.GT.1) THEN
          DO iptfr=1,nptir
            i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
!           AVAILABLE DEPTH AND TRACER QUANTITY ARE SHARED BETWEEN PROCESSORS
!           FOR POINTS GIVING WATER, THEY WILL BE CANCELLED LATER
            ht%R(i)=ht%R(i)*mesh%IFAC%I(i)
            t5%R(i)=ht%R(i)*f%R(i)
          ENDDO
        ENDIF
      ELSE
!       NOT ALL THE POINTS NEED TO BE INITIALISED NOW
        DO ir=1,remain
          i=indic_cpos(ir)
          i1=gloseg1(i)
          i2=gloseg2(i)
          t1%R(i1)=0.d0
          t1%R(i2)=0.d0
!         SAVING THE DEPTH AND TRACER
          t4%R(i1)=ht%R(i1)
          t4%R(i2)=ht%R(i2)
          t6%R(i1)=f%R(i1)
          t6%R(i2)=f%R(i2)
          t5%R(i1)=ht%R(i1)*f%R(i1)
          t5%R(i2)=ht%R(i2)*f%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 AND TRACER
            t4%R(i)=ht%R(i)
            t6%R(i)=f%R(i)
!           AVAILABLE DEPTH AND TRACER QUANTITY ARE SHARED BETWEEN PROCESSORS
!           FOR POINTS GIVING WATER, THEY WILL BE CANCELLED LATER
            ht%R(i)=ht%R(i)*mesh%IFAC%I(i)
            t5%R(i)=ht%R(i)*f%R(i)
          ENDDO
        ENDIF
      ENDIF
      DO ir=1,remain
        i=indic_cpos(ir)
        i1=gloseg1(i)
        i2=gloseg2(i)
        IF(fxmat(i).GT.eps_flux) THEN
          t1%R(i1)=t1%R(i1)+fxmat(i)
          ht%R(i1)=0.d0
          t5%R(i1)=0.d0
        ELSEIF(fxmat(i).LT.-eps_flux) THEN
          t1%R(i2)=t1%R(i2)-fxmat(i)
          ht%R(i2)=0.d0
          t5%R(i2)=0.d0
        ENDIF
      ENDDO
!
      IF(ncsize.GT.1) CALL parcom(t1,2,mesh)
!
!     FOR ISOLATED POINTS CONNECTED TO AN ACTIVE SEGMENT
!     THAT IS IN ANOTHER SUBDOMAIN
      IF(ncsize.GT.1) THEN
        DO iptfr=1,nptir
          i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
          IF(t1%R(i).GT.eps_flux) THEN
            ht%R(i)=0.d0
            t5%R(i)=0.d0
          ENDIF
        ENDDO
      ENDIF
!
      c=0.d0
      newremain=0
!
      DO ir=1,remain
        i=indic_cpos(ir)
        i1=gloseg1(i)
        i2=gloseg2(i)
        IF(fxmat(i).GT.eps_flux) THEN
!         SHARING ON DEMAND: FRACTION OF DEPTH TAKEN
!         T4 IS THE STORED DEPTH
!         1.D-20 ADDED HERE OTHERWISE IT OCCURED THAT THE PRODUCT OF 2
!         STRICTLY POSITIVE NUMBERS GAVE 0.D0
          IF(t4%R(i1).GT.1.d-20) THEN
            hseg1=t4%R(i1)*fxmat(i)/t1%R(i1)
!           END OF SHARING ON DEMAND
            hfl1= dt*unsv2d%R(i1)*fxmat(i)
            IF(hfl1.GT.hseg1) THEN
              tet=hseg1/hfl1
!             HSEG2 AND THUS HT WILL BE STRICTLY POSITIVE
              hseg2=dt*unsv2d%R(i2)*fxmat(i)*tet
              ht%R(i2)=ht%R(i2)+hseg2
!             GROUPING H*F
              t5%R(i2)=t5%R(i2)+hseg2*t6%R(i1)
!             RECOMPUTING F (AS WEIGHTED AVERAGE)
!             THIS MAY BE DONE SEVERAL TIMES FOR THE SAME POINT
!             BUT THE LAST ONE WILL BE THE GOOD ONE
              f%R(i2)=t5%R(i2)/ht%R(i2)
              fxmat(i)=fxmat(i)*(1.d0-tet)
              c=c+fxmat(i)
              newremain=newremain+1
              indic_cpos(newremain)=i
            ELSE
              hseg1=hseg1-hfl1
              hseg2=dt*unsv2d%R(i2)*fxmat(i)
              ht%R(i2)=ht%R(i2)+hseg2
              t5%R(i2)=t5%R(i2)+hseg2*t6%R(i1)
!             THE LAST ONE WILL BE THE GOOD ONE
              f%R(i2)=t5%R(i2)/ht%R(i2)
              IF(hseg1.GT.0.d0) THEN
                ht%R(i1)=ht%R(i1)+hseg1
                t5%R(i1)=t5%R(i1)+hseg1*t6%R(i1)
!               THE LAST ONE WILL BE THE GOOD ONE
                f%R(i1)=t5%R(i1)/ht%R(i1)
              ENDIF
            ENDIF
          ELSE
!           NO WATER NO FLUX TRANSMITTED, NOTHING CHANGED
            c=c+fxmat(i)
            newremain=newremain+1
            indic_cpos(newremain)=i
          ENDIF
        ELSEIF(fxmat(i).LT.-eps_flux) THEN
!         SHARING ON DEMAND
          IF(t4%R(i2).GT.1.d-20) THEN
!         1.D-20 ADDED HERE OTHERWISE IT OCCURED THAT THE PRODUCT OF 2
!         STRICTLY POSITIVE NUMBERS GAVE 0.D0
            hseg2=-t4%R(i2)*fxmat(i)/t1%R(i2)
!           END OF SHARING ON DEMAND
            hfl2=-dt*unsv2d%R(i2)*fxmat(i)
            IF(hfl2.GT.hseg2) THEN
              tet=hseg2/hfl2
!             HSEG1 AND THUS HT WILL BE STRICTLY POSITIVE
              hseg1=-dt*unsv2d%R(i1)*fxmat(i)*tet
              ht%R(i1)=ht%R(i1)+hseg1
              t5%R(i1)=t5%R(i1)+hseg1*t6%R(i2)
!             THE LAST ONE WILL BE THE GOOD ONE
              f%R(i1)=t5%R(i1)/ht%R(i1)
              fxmat(i)=fxmat(i)*(1.d0-tet)
              c=c-fxmat(i)
              newremain=newremain+1
              indic_cpos(newremain)=i
            ELSE
              hseg1=-dt*unsv2d%R(i1)*fxmat(i)
              hseg2=hseg2-hfl2
              ht%R(i1)=ht%R(i1)+hseg1
              t5%R(i1)=t5%R(i1)+hseg1*t6%R(i2)
              f%R(i1)=t5%R(i1)/ht%R(i1)
              IF(hseg2.GT.0.d0) THEN
                ht%R(i2)=ht%R(i2)+hseg2
                t5%R(i2)=t5%R(i2)+hseg2*t6%R(i2)
!               THE LAST ONE WILL BE THE GOOD ONE
                f%R(i2)=t5%R(i2)/ht%R(i2)
              ENDIF
            ENDIF
          ELSE
!           NO WATER NO FLUX TRANSMITTED, NOTHING CHANGED
            c=c-fxmat(i)
            newremain=newremain+1
            indic_cpos(newremain)=i
          ENDIF
        ENDIF
      ENDDO
!
!     MERGING DEPTHS AND F AT INTERFACE POINTS
!
      IF(ncsize.GT.1) THEN
        DO iptfr=1,nptir
!         ARRAY WITH HT*F AT INTERFACE POINTS
          i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
          t1%R(i)=ht%R(i)*f%R(i)
        ENDDO
!       SUMMING HT*F AT INTERFACE POINTS
        CALL parcom(t1,2,mesh)
!       SUMMING THE NEW POSITIVE PARTIAL DEPTHS OF INTERFACE POINTS
        CALL parcom(ht,2,mesh)
!       AVERAGE F AT INTERFACE POINTS
        DO iptfr=1,nptir
          i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
          IF(ht%R(i).GT.0.d0) f%R(i)=t1%R(i)/ht%R(i)
        ENDDO
      ENDIF
!
      IF(ncsize.GT.1) c=p_sum(c)
      IF(testing) WRITE(lu,*) 'FLUX NON PRIS EN COMPTE=',c
!
      remain=newremain
!
      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 NEGATIVE SOURCES
!
      IF(yasmh) THEN
        IF(optsou.EQ.1) THEN
          DO i=1,npoin
            IF(smh%R(i).LT.0.d0) THEN
              ht%R(i)=ht%R(i)+dt*smh%R(i)
              f%R(i)=f%R(i)+dt/max(ht%R(i),1.d-4)*smh%R(i)*fscexp%R(i)
            ENDIF
          ENDDO
        ELSEIF(optsou.EQ.2) THEN
          DO i=1,npoin
            IF(smh%R(i).LT.0.d0) THEN
              ht%R(i)=ht%R(i)+dt*smh%R(i)*unsv2d%R(i)
              f%R(i)=f%R(i)+dt/max(ht%R(i),1.d-4)*
     &                        unsv2d%R(i)*smh%R(i)*fscexp%R(i)
            ENDIF
          ENDDO
        ENDIF
      ENDIF
!
!     RAIN-EVAPORATION: RAIN DONE ABOVE, NOW EVAPORATION
!
      IF(rain) THEN
        DO i=1,npoin
          c=min(pluie%R(i),0.d0)
!         POSITIVITY NOT TESTED HERE, WOULD REQUIRE C=MAX(C,-HT%R(I)/DT)
!         BUT THEN MASS-BALANCE WOULD NOT BE CORRECT,
          ht%R(i)=ht%R(i)+dt*c
!                                                VALUE IN VAPOR
          f%R(i)=f%R(i)+dt/max(ht%R(i),1.d-4)*c*(0.d0-f%R(i))
        ENDDO
      ENDIF
!
!     BOUNDARY FLUXES : ADDING THE EXITING (POSITIVE) FLUXES
!                       WITH A POSSIBLE LIMITATION
!
      DO iptfr=1,nptfr
        i=nbor(iptfr)
!                               T2 = // ASSEMBLED FLBOR
        hfl1=dt*unsv2d%R(i)*max(t2%R(i),0.d0)
        tet=1.d0
        IF(hfl1.GT.ht%R(i)) tet=ht%R(i)/hfl1
!       MAX IS ONLY TO PREVENT TRUNCATION ERROR
        ht%R(i)=max(ht%R(i)-hfl1*tet,0.d0)
!       LIMITATION OF FLBOR (MUST HAVE BEEN DONE ALREADY
!                            IN POSITIVE_DEPTHS)
!       FLBOR%R(IPTFR)=FLBOR%R(IPTFR)*TET
        IF(limtra(iptfr).EQ.kdir) THEN
          f%R(i)=f%R(i)-hfl1*tet/max(ht%R(i),1.d-4)*
     &           (fbor%R(iptfr)-f%R(i))
          flbortra%R(iptfr)=flbor%R(iptfr)*fbor%R(iptfr)
        ELSEIF(limtra(iptfr).EQ.kddl) THEN
          flbortra%R(iptfr)=flbor%R(iptfr)*f%R(i)
        ENDIF
      ENDDO
!
      IF(testing) THEN
        c=0.d0
        DO i=1,npoin
          c=c+(ht%R(i)-h%R(i))**2
        ENDDO
!                       FAUX MAIS PAS GRAVE SI 0.
        IF(ncsize.GT.1) c=p_sum(c)
        WRITE(lu,*) 'DIFFERENCE ENTRE H ET HT =',c
!
        c=1.d99
        DO i=1,npoin
          c=min(c,f%R(i))
        ENDDO
        IF(ncsize.GT.1) c=p_min(c)
        WRITE(lu,*) 'APRES TRAITEMENT TRACEUR MIN=',c
        c=-1.d99
        DO i=1,npoin
          c=max(c,f%R(i))
        ENDDO
        IF(ncsize.GT.1) c=p_max(c)
        WRITE(lu,*) 'APRES TRAITEMENT TRACEUR MAX=',c
      ENDIF
!
!-----------------------------------------------------------------------
!
!     SOURCE TERMS (TREATED OUTSIDE THE SUB-ITERATION LOOP, HENCE
!                   ORIGINAL TIME-STEP DT0)
!
      IF(yasmi) THEN
!
!       IMPLICIT AND EXPLICIT SOURCE TERM
!
        IF(iopt2.EQ.0) THEN
          DO i = 1,mesh%NPOIN
            f%R(i)=(f%R(i)+dt*sm%R(i))/
     &           (1.d0-dt*smi%R(i)/max(h%R(i),1.d-15))
!           COULD BE DONE LIKE THIS...
!           F%R(I)=H%R(I)*(F%R(I)+DT*SM%R(I))/(H%R(I)-DT0*SMI%R(I))
          ENDDO
        ELSEIF(iopt2.EQ.1) THEN
!         HERE WE ASSUME THAT SMI WILL PREVENT A DIVISION BY ZERO
!         THIS IS THE CASE WITH SETTLING VELOCITY IN SISYPHE
          DO i = 1,mesh%NPOIN
          f%R(i)=(f%R(i)*ht%R(i)+dt*sm%R(i)*h%R(i))/
     &           (h%R(i)-dt*smi%R(i))
          ENDDO
        ENDIF
!
      ELSE
!
!       EXPLICIT SOURCE TERM ONLY (AND IOPT2=1 NOT TREATED !!!)
!
        DO i = 1,mesh%NPOIN
          f%R(i) = f%R(i)+dt*sm%R(i)
        ENDDO
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(entet) THEN
        WRITE(lu,202) niter
      ENDIF
!
      IF(niter.EQ.nitmax) THEN
        WRITE(lu,203) niter
      ENDIF
!
202   FORMAT('CVTRVF_NERD (SCHEME NERD, 13 OR 14): ',1i3,' ITERATIONS')
203   FORMAT(' CVTRVF_NERD (SCHEME NERD, 13 OR 14): ',1i3,' ITERATIONS',
     &  ' = MAXIMUM',
     &  /,'INCREASE MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES')
!
!-----------------------------------------------------------------------
!
      RETURN
      END