cvtrvf_nerd_2.f

Go to the documentation of this file.

!                   ************************
                    SUBROUTINE cvtrvf_nerd_2
!                   ************************
!
     &(f1,f1n,f1scexp,f2,f2n,f2scexp,h,hn,hprop,udel,vdel,dm1,
     & zconv,solsys,sm1,sm2,smh,yasmh,smi1,smi2,yasmi,
     & f1bor,f2bor,masktr,mesh,t1,t2,t3,t4,t5,t6,t7,t8,ht,
     & dt,entet,msk,maskel,optsou,
     & limtra1,limtra2,kdir,kddl,nptfr,flbor,yaflbor,unsv2d,iopt,
     & flbortra1,flbortra2,gloseg1,gloseg2,nbor,
     & rain,pluie,train1,train2,nitmax)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    FINITE VOLUMES, UPWIND, EXPLICIT AND MONOTONIC
!+               ADVECTOR EVEN WITH TIDAL FLATS.
!+               THIS IS A COPY OF CVTRVF_POS, WRITTEN FOR 2 VARIABLES.
!
!warning  AFBOR AND BFBOR MUST BE 0 FOR THE BOUNDARY ELEMENTS
!+            WITH NO FRICTION
!warning  DISCRETISATION OF VISC
!
!history  J-M HERVOUET   (LNHE)
!+        19/11/2010
!+        V6P0
!+   OPTIMIZATION (2 ABS SUPPRESSED)
!
!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  J-M HERVOUET (LNHE)
!+        12/07/2012
!+        V6P2
!+   T2 set to 0 to start with required in parallel
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        16/04/2013
!+        V6P2
!+   Intent of FLBOR changed and conditional building of FLBOR added.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        30/05/2013
!+        V6P2
!+   Argument NITMAX added.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        12/06/2015
!+        V7P1
!+   Adaptation to the fact that MESH%FAC is now replaced by MESH%IFAC.
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| 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.
!| F1             |<--| F1 AT TIME T(N+1)
!| F2             |<--| F2 AT TIME T(N+1)
!| F1BOR          |-->| DIRICHLET CONDITIONS ON F1.
!| F2BOR          |-->| DIRICHLET CONDITIONS ON F2.
!| FLBOR          |-->| FLUXES AT BOUNDARIES
!| FLBORTRA       |<->| TRACER FLUXES AT BOUNDARIES
!| F1N            |-->| F1 AT TIME T(N)
!| F2N            |-->| F2 AT TIME T(N)
!| F1SCEXP        |-->| EXPLICIT PART OF THE F1 SOURCE TERM
!|                |   | EQUAL TO ZERO EVERYWHERE BUT ON SOURCES
!|                |   | WHERE THERE IS FSCE - (1-TETAT) FN
!|                |   | SEE DIFSOU
!| F2SCEXP        |-->| EXPLICIT PART OF THE F1 SOURCE TERM
!| GLOSEG1        |-->| FIRST POINT OF SEGMENTS
!| GLOSEG2        |-->| SECOND POINT OF SEGMENTS
!| HT             |<--| WORK ARRAYS (MODIFIED DEPTHS TO TAKE MASS-LUMPING
!|                |   | INTO ACCOUNT)
!| 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
!| LIMTRA1        |-->| BOUNDARY CONDITIONS OF F1 ON BOUNDARY POINTS
!| LIMTRA2        |-->| BOUNDARY CONDITIONS OF F2 ON BOUNDARY 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
!| PLUIE          |-->| RAIN OR EVAPORATION IN MM/S IN A BIEF_OBJ
!| RAIN           |-->| IF YES, RAIN OR EVAPORATION
!| SM1            |-->| SOURCE TERMS OF F1.
!| SM2            |-->| SOURCE TERMS OF F2.
!| SMH            |-->| SOURCE TERM IN CONTINUITY EQUATION
!| SMI1           |-->| IMPLICIT SOURCE TERM OF F1.
!| SMI2           |-->| IMPLICIT SOURCE TERM OF F2.
!| 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
!| T8             |<->| WORK BIEF_OBJ STRUCTURE
!| TRAIN1         |-->| VALUE OF TRACER 1 IN THE RAIN
!| TRAIN2         |-->| VALUE OF TRACER 2 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_2 => cvtrvf_nerd_2
      USE declarations_telemac, ONLY : deja_cpos2, indic_cpos2
!
      USE declarations_special
      USE interface_parallel, ONLY : p_sum,p_min,p_max
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: OPTSOU,KDIR,NPTFR,SOLSYS
      INTEGER, INTENT(IN)             :: KDDL,IOPT,NITMAX
      INTEGER, INTENT(IN)             :: GLOSEG1(*),GLOSEG2(*)
      INTEGER, INTENT(IN)             :: LIMTRA1(nptfr),NBOR(nptfr)
      INTEGER, INTENT(IN)             :: LIMTRA2(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DT
      DOUBLE PRECISION, INTENT(IN)    :: TRAIN1,TRAIN2
      LOGICAL, INTENT(IN)             :: YASMH,YAFLBOR,RAIN
      LOGICAL, INTENT(IN)             :: MSK,ENTET,YASMI
      TYPE(bief_obj), INTENT(IN)      :: MASKEL,H,HN,DM1,ZCONV
      TYPE(bief_obj), INTENT(IN)      :: UNSV2D,HPROP
      TYPE(bief_obj), INTENT(INOUT)   :: F1,SM1,F2,SM2,HT
      TYPE(bief_obj), INTENT(IN)      :: F1BOR,UDEL,VDEL,F1N,SMI1,SMH
      TYPE(bief_obj), INTENT(IN)      :: F2BOR,F2N,SMI2
      TYPE(bief_obj), INTENT(INOUT)   :: FLBORTRA1,FLBORTRA2
      TYPE(bief_obj), INTENT(INOUT)   :: T1,T2,T3,T4,T5,T6,T7,T8
      TYPE(bief_obj), INTENT(IN)      :: F1SCEXP,F2SCEXP,MASKTR
      TYPE(bief_obj), INTENT(INOUT)   :: FLBOR
      TYPE(bief_obj), INTENT(IN)      :: PLUIE
      TYPE(bief_mesh)                 :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,IOPT1,IOPT2,NPOIN,IPTFR,I1,I2,NITER,REMAIN_SEG,NEWREMAIN
      INTEGER IR
!
!-----------------------------------------------------------------------
!
      DOUBLE PRECISION C,CPREV,CINIT,HFL1,HFL2,TET
      DOUBLE PRECISION HSEG1,HSEG2
      CHARACTER(LEN=16) FORMUL
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: FXMAT
      LOGICAL, PARAMETER :: TESTING = .false.
      DOUBLE PRECISION, PARAMETER :: EPS_FLUX = 1.d-15
!
!-----------------------------------------------------------------------
!
!     INDIC_CPOS2 WILL BE A LIST OF SEGMENTS WITH NON ZERO FLUXES
!
      IF(.NOT.deja_cpos2) THEN
        ALLOCATE(indic_cpos2(mesh%NSEG))
        deja_cpos2=.true.
      ENDIF
!
!-----------------------------------------------------------------------
!
      fxmat=>mesh%MSEG%X%R(1:mesh%NSEG)
!
!-----------------------------------------------------------------------
!
      npoin=h%DIM1
!
!     EXTRACTING OPTIONS, AND CONTROL
!
      iopt2=iopt/10
      iopt1=iopt-10*iopt2
      IF(iopt1.LT.0.OR.iopt1.GT.3) THEN
        WRITE(lu,*) 'CVTRVF_NERD_2: OPTION IOPT1 UNKNOWN: ',iopt1
        CALL plante(1)
        stop
      ENDIF
      IF(iopt2.NE.0.AND.iopt2.NE.1) THEN
        WRITE(lu,*) 'CVTRVF_NERD_2: 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
!
!     CALCUL DES FLUX PAR NOEUDS
!
      formul='HUGRADP         '
      IF(solsys.EQ.2) formul(8:8)='2'
      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
!     CALCUL DES FLUX PAR SEGMENT (TE1 SUIVI DE FALSE NON UTILISE)
!     FXMAT IS NOT ASSEMBLED IN //
!
!----------------------------------------
! DIFFERENT OPTIONS TO COMPUTE THE FLUXES
!----------------------------------------
!
      CALL flux_ef_vf(fxmat,mesh%W%R,mesh%NSEG,mesh%NELEM,mesh%NELMAX,
     &                mesh%ELTSEG%I,mesh%ORISEG%I,
     &                mesh%IKLE%I,.true.,iopt1)
!
!----------------------------------------
!
!     AVERAGING FLUXES ON INTERFACE SEGMENTS BY ASSEMBLING AND
!     DIVIDING BY 2. THIS WILL GIVE THE UPWINDING INFORMATION
!
      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
!
!----------------------------------------
! END OF THE OPTIONS
!----------------------------------------
!
      CALL cpstvc(h,t2)
!     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)
!
!     INITIALIZING F1 AND F2 AT THE OLD VALUE
!
      CALL os('X=Y     ',x=f1,y=f1n)
      CALL os('X=Y     ',x=f2,y=f2n)
!
      cprev=0.d0
      DO i=1,mesh%NSEG
        cprev=cprev+abs(fxmat(i))
      ENDDO
      IF(ncsize.GT.1) cprev=p_sum(cprev)
      cinit=cprev
      IF(testing) WRITE(lu,*) 'SOMME INITIALE DES FLUX=',cprev
!
!     BOUCLE SUR LES SEGMENTS, POUR PRENDRE EN COMPTE LES FLUX
!     ADMISSIBLES
!
!     ADDING THE SOURCES (SMH IS NATURALLY ASSEMBLED IN //)
      IF(yasmh) THEN
        IF(optsou.EQ.1) THEN
          DO i=1,npoin
            ht%R(i)=hn%R(i)+dt*smh%R(i)
            f1%R(i)=f1n%R(i)+dt/max(ht%R(i),1.d-4)*
     &                       smh%R(i)*f1scexp%R(i)
            f2%R(i)=f2n%R(i)+dt/max(ht%R(i),1.d-4)*
     &                       smh%R(i)*f2scexp%R(i)
          ENDDO
        ELSEIF(optsou.EQ.2) THEN
          DO i=1,npoin
            ht%R(i)=hn%R(i)+dt*smh%R(i)*unsv2d%R(i)
            f1%R(i)=f1n%R(i)+dt/max(ht%R(i),1.d-4)*
     &                       unsv2d%R(i)*smh%R(i)*f1scexp%R(i)
            f2%R(i)=f2n%R(i)+dt/max(ht%R(i),1.d-4)*
     &                       unsv2d%R(i)*smh%R(i)*f2scexp%R(i)
          ENDDO
        ENDIF
      ELSE
        DO i=1,npoin
          ht%R(i)=hn%R(i)
        ENDDO
      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
          f1%R(i)=f1%R(i)+dt/max(ht%R(i),1.d-4)*c*(train1-f1%R(i))
          f2%R(i)=f2%R(i)+dt/max(ht%R(i),1.d-4)*c*(train2-f2%R(i))
        ENDDO
      ENDIF
!
      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)
      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
!       THE FINAL DEPTH IS TAKEN
        IF(limtra1(iptfr).EQ.kdir) THEN
          f1%R(i)=f1n%R(i)-dt/max(ht%R(i),1.d-4)*
     &    unsv2d%R(i)*min(t2%R(i),0.d0)*(f1bor%R(iptfr)-f1n%R(i))
        ELSEIF(limtra1(iptfr).EQ.kddl) THEN
          IF(t2%R(i).LE.0.d0) THEN
!           FLBORTRA1 IS NOT ASSEMBLED
            flbortra1%R(iptfr)=flbor%R(iptfr)*f1n%R(i)
          ENDIF
        ENDIF
        IF(limtra2(iptfr).EQ.kdir) THEN
          f2%R(i)=f2n%R(i)-dt/max(ht%R(i),1.d-4)*
     &    unsv2d%R(i)*min(t2%R(i),0.d0)*(f2bor%R(iptfr)-f2n%R(i))
        ELSEIF(limtra2(iptfr).EQ.kddl) THEN
          IF(t2%R(i).LE.0.d0) THEN
!           FLBORTRA2 IS NOT ASSEMBLED
            flbortra2%R(iptfr)=flbor%R(iptfr)*f1n%R(i)
          ENDIF
        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_seg=mesh%NSEG
      DO i=1,remain_seg
        indic_cpos2(i)=i
      ENDDO
!
      niter = 0
777   CONTINUE
      niter = niter + 1
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     FOR DISTRIBUTING THE DEPTHS BETWEEN SEGMENTS
!
!     T1 : TOTAL FLUX REMOVED OF EACH POINT
!     T4 : DEPTH H SAVED
!     T5 : H*F1 SAVED
!     T6 : F1 SAVED
!     T7 : H*F2 SAVED
!     T8 : F2 SAVED
!
      CALL cpstvc(h,t1)
      IF(niter.EQ.1) THEN
        DO i=1,npoin
          t1%R(i)=0.d0
          t4%R(i)=ht%R(i)
          t6%R(i)=f1%R(i)
          t8%R(i)=f2%R(i)
          t5%R(i)=ht%R(i)*f1%R(i)
          t7%R(i)=ht%R(i)*f2%R(i)
        ENDDO
        IF(ncsize.GT.1) THEN
          DO iptfr=1,nptir
            i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
!           AVAILABLE DEPTH IS SHARED BETWEEN PROCESSORS
            ht%R(i)=ht%R(i)*mesh%IFAC%I(i)
            t5%R(i)=ht%R(i)*f1%R(i)
            t7%R(i)=ht%R(i)*f2%R(i)
          ENDDO
        ENDIF
      ELSE
!       NOT ALL THE POINTS NEED TO BE INITIALISED NOW
        DO ir=1,remain_seg
          i=indic_cpos2(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)=f1%R(i1)
          t6%R(i2)=f1%R(i2)
          t8%R(i1)=f2%R(i1)
          t8%R(i2)=f2%R(i2)
          t5%R(i1)=ht%R(i1)*f1%R(i1)
          t5%R(i2)=ht%R(i2)*f1%R(i2)
          t7%R(i1)=ht%R(i1)*f2%R(i1)
          t7%R(i2)=ht%R(i2)*f2%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)=f1%R(i)
            t8%R(i)=f2%R(i)
!           AVAILABLE DEPTH IS SHARED BETWEEN PROCESSORS
            ht%R(i)=ht%R(i)*mesh%IFAC%I(i)
            t5%R(i)=ht%R(i)*f1%R(i)
            t7%R(i)=ht%R(i)*f2%R(i)
          ENDDO
        ENDIF
      ENDIF
      DO ir=1,remain_seg
        i=indic_cpos2(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
          t7%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
          t7%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
            t7%R(i)=0.d0
          ENDIF
        ENDDO
      ENDIF
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      c=0.d0
      newremain=0
!
      DO ir=1,remain_seg
        i=indic_cpos2(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
          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)
              t7%R(i2)=t7%R(i2)+hseg2*t8%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
              f1%R(i2)=t5%R(i2)/ht%R(i2)
              f2%R(i2)=t7%R(i2)/ht%R(i2)
              fxmat(i)=fxmat(i)*(1.d0-tet)
              c=c+fxmat(i)
              newremain=newremain+1
              indic_cpos2(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)
              t7%R(i2)=t7%R(i2)+hseg2*t8%R(i1)
!             THE LAST ONE WILL BE THE GOOD ONE
              f1%R(i2)=t5%R(i2)/ht%R(i2)
              f2%R(i2)=t7%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)
                t7%R(i1)=t7%R(i1)+hseg1*t8%R(i1)
!               THE LAST ONE WILL BE THE GOOD ONE
                f1%R(i1)=t5%R(i1)/ht%R(i1)
                f2%R(i1)=t7%R(i1)/ht%R(i1)
              ENDIF
            ENDIF
          ELSE
!           NO WATER NO FLUX TRANSMITTED, NOTHING CHANGED
            c=c+fxmat(i)
            newremain=newremain+1
            indic_cpos2(newremain)=i
          ENDIF
        ELSEIF(fxmat(i).LT.-eps_flux) THEN
!         SHARING ON DEMAND
          IF(t4%R(i2).GT.1.d-20) THEN
            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)
              t7%R(i1)=t7%R(i1)+hseg1*t8%R(i2)
!             THE LAST ONE WILL BE THE GOOD ONE
              f1%R(i1)=t5%R(i1)/ht%R(i1)
              f2%R(i1)=t7%R(i1)/ht%R(i1)
              fxmat(i)=fxmat(i)*(1.d0-tet)
              c=c-fxmat(i)
              newremain=newremain+1
              indic_cpos2(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)
              t7%R(i1)=t7%R(i1)+hseg1*t8%R(i2)
              f1%R(i1)=t5%R(i1)/ht%R(i1)
              f2%R(i1)=t7%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)
                t7%R(i2)=t7%R(i2)+hseg2*t8%R(i2)
!               THE LAST ONE WILL BE THE GOOD ONE
                f1%R(i2)=t5%R(i2)/ht%R(i2)
                f2%R(i2)=t7%R(i2)/ht%R(i2)
              ENDIF
            ENDIF
          ELSE
!           NO WATER NO FLUX TRANSMITTED, NOTHING CHANGED
            c=c-fxmat(i)
            newremain=newremain+1
            indic_cpos2(newremain)=i
          ENDIF
        ENDIF
      ENDDO
!
      remain_seg=newremain
!
!     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)*f1%R(i)
          t3%R(i)=ht%R(i)*f2%R(i)
        ENDDO
!       SUMMING HT*F AT INTERFACE POINTS
        CALL parcom(t1,2,mesh)
        CALL parcom(t3,2,mesh)
!       SUMMING THE NEW POSITIVE PARTIAL DEPTHS OF INTERFACE POINTS
        CALL parcom(ht,2,mesh)
!       AVERAGE F1 AND F2 AT INTERFACE POINTS
        DO iptfr=1,nptir
          i=mesh%NACHB%I(nbmaxnshare*(iptfr-1)+1)
          IF(ht%R(i).GT.0.d0) THEN
            f1%R(i)=t1%R(i)/ht%R(i)
            f2%R(i)=t3%R(i)/ht%R(i)
          ENDIF
        ENDDO
      ENDIF
!
      IF(ncsize.GT.1) c=p_sum(c)
      IF(testing) WRITE(lu,*) 'FLUX NON PRIS EN COMPTE=',c
      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
!
!     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
!         F1%R(I)=F1%R(I)+DT/MAX(HT%R(I),1.D-4)*C*(0.D0-F1%R(I))
          f1%R(i)=f1%R(i)-dt/max(ht%R(i),1.d-4)*c*f1%R(i)
!         F2%R(I)=F2%R(I)+DT/MAX(HT%R(I),1.D-4)*C*(0.D0-F2%R(I))
          f2%R(i)=f2%R(i)-dt/max(ht%R(i),1.d-4)*c*f2%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
!       NEXT LINE SHOULD NEVER HAPPEN (DONE IN POSITIVE_DEPTHS)
        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(limtra1(iptfr).EQ.kdir) THEN
          f1%R(i)=f1%R(i)-hfl1*tet/max(ht%R(i),1.d-4)*
     &           (f1bor%R(iptfr)-f1%R(i))
          flbortra1%R(iptfr)=flbor%R(iptfr)*f1bor%R(iptfr)
        ELSEIF(limtra1(iptfr).EQ.kddl) THEN
          IF(t2%R(i).GT.0.d0) THEN
            flbortra1%R(iptfr)=flbor%R(iptfr)*f1%R(i)
          ENDIF
        ELSE
          flbortra1%R(iptfr)=0.d0
        ENDIF
        IF(limtra2(iptfr).EQ.kdir) THEN
          f2%R(i)=f2%R(i)-hfl1*tet/max(ht%R(i),1.d-4)*
     &           (f2bor%R(iptfr)-f2%R(i))
          flbortra2%R(iptfr)=flbor%R(iptfr)*f2bor%R(iptfr)
        ELSEIF(limtra2(iptfr).EQ.kddl) THEN
          IF(t2%R(i).GT.0.d0) THEN
            flbortra2%R(iptfr)=flbor%R(iptfr)*f2%R(i)
          ENDIF
        ELSE
          flbortra2%R(iptfr)=0.d0
        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,f1%R(i))
        ENDDO
        IF(ncsize.GT.1) c=p_min(c)
        WRITE(lu,*) 'APRES TRAITEMENT TRACEUR 1 MIN=',c
        c=-1.d99
        DO i=1,npoin
          c=max(c,f1%R(i))
        ENDDO
        IF(ncsize.GT.1) c=p_max(c)
        WRITE(lu,*) 'APRES TRAITEMENT TRACEUR 1 MAX=',c
      ENDIF
!
!-----------------------------------------------------------------------
!
!     EXPLICIT SOURCE TERM
!
      DO i = 1,mesh%NPOIN
        f1%R(i) = f1%R(i)+dt*sm1%R(i)
        f2%R(i) = f2%R(i)+dt*sm2%R(i)
      ENDDO
!
!     IMPLICIT SOURCE TERM
!
      IF(yasmi) THEN
        DO i = 1,mesh%NPOIN
          f1%R(i) = f1%R(i)/(1.d0-dt*smi1%R(i)/max(h%R(i),1.d-15))
          f2%R(i) = f2%R(i)/(1.d0-dt*smi2%R(i)/max(h%R(i),1.d-15))
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(entet) THEN
        WRITE(lu,102) niter
      ENDIF
!
      IF(niter.EQ.nitmax) THEN
        WRITE(lu,103) niter
      ENDIF
!
102   FORMAT(' CVTRVF_NERD_2 (SCHEME 13 OR 14): ',1i3,' ITERATIONS')
103   FORMAT(' CVTRVF_NERD_2 (SCHEME NERD, 13 OR 14): ',1i3,
     &  ' ITERATIONS = MAXIMUM',
     &  /,'INCREASE MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES')
!
!-----------------------------------------------------------------------
!
      RETURN
      END