cvtrvf_eria.f

Go to the documentation of this file.

!                   **********************
                    SUBROUTINE cvtrvf_eria
!                   **********************
!
     &(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,nbor,
     & rain,pluie,train,nitmax,nco_dist,optadv)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    Distributive advection scheme ERIA
!
!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   (EDF LAB, LNHE)
!+        09/09/2016
!+        V7P2
!+   First version. Comes from the former cvtrvf_pos with option 1.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        23/11/2016
!+        V7P3
!+   Adding an option for predictor for sharing volumes of points
!+   between elements (see variable OPTPRE, OPTPRE=1 old strategy,
!+   OPTPRE=2, new and simpler strategy). Tests do not show much
!+   difference, kept here to be tested in other cases.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        02/01/2017
!+        V7P3
!+   Second order in time added.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        09/09/2017
!+        V7P3
!+   Treating cases where NELMAX is not equal to NELEM.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| 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
!| 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
!| 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
!| NCO_DIST       |-->| NUMBER OF CORRECTIONS IN DISTRIBUTIVE SCHEMES
!| NITMAX         |-->| MAXIMUM NUMBER OF ITERATIONS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| OPTADV         |-->| FOR ERIA SCHEME 1: FIRST ORDER
!|                |   |                 2: SECOND ORDER
!| 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_eria => cvtrvf_eria
      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,OPTADV
      INTEGER, INTENT(IN)             :: NCO_DIST
      INTEGER, INTENT(IN)             :: NBOR(nptfr)
      INTEGER, INTENT(INOUT)          :: LIMTRA(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DT,TRAIN
      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(IN)      :: UDEL,VDEL,FN,SMI,SMH
      TYPE(bief_obj), INTENT(INOUT)   :: FLBORTRA,FBOR,F,SM,HT
      TYPE(bief_obj), INTENT(INOUT)   :: T1,T2,T3,T4,T5,T6,T7,FLBOR
      TYPE(bief_obj), INTENT(IN)      :: FSCEXP,MASKTR,PLUIE
      TYPE(bief_mesh)                 :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,IOPT2,NPOIN,IPTFR,I1,I2,I3,NITER,NEWREMAIN,NELMAX
      INTEGER IR,N,NELEM,IELEM,REMAIN,NSEG
      INTEGER OPTPRE
!
!-----------------------------------------------------------------------
!
      DOUBLE PRECISION C,CPREV,CINIT,HFL1,TET,LOCALMIN,LOCALMAX
      DOUBLE PRECISION F1,F2,F3,COEF,FI1,FI2,FI3,VOL1,VOL2,VOL3
      DOUBLE PRECISION DT1,DT2,DT3,FMIN,FMAX,VOLD1,VOLD2,VOLD3
      DOUBLE PRECISION SURDT,FITOT,BETA1,BETA2,BETA3,A1,A2,A3
      DOUBLE PRECISION FP1,FP2,FP3,FIP1,FIP2,FIP3,FS1,FS2,FS3
      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
      DOUBLE PRECISION, POINTER :: FLOP1(:),FLOP2(:),FLOP3(:)
      DOUBLE PRECISION, POINTER :: DTLIM1(:),DTLIM2(:),DTLIM3(:)
      DOUBLE PRECISION, POINTER :: SVOL1(:),SVOL2(:),SVOL3(:)
!
!-----------------------------------------------------------------------
!
!     INDIC_CPOS WILL BE A LIST OF SEGMENTS WITH NON ZERO FLUXES
!
      nelem=mesh%NELEM
      nelmax=mesh%NELMAX
      nseg=mesh%NSEG
      npoin=h%DIM1
      surdt=1.d0/dt
!
      IF(.NOT.deja_cpos) THEN
        ALLOCATE(indic_cpos(nelem))
        deja_cpos=.true.
      ENDIF
!
!-----------------------------------------------------------------------
!
!     THERE IS PLENTY OF MEMORY AVAILABLE IN A BIEF_MESH STRUCTURE...
!
!     3*NELEM+NPTFR=2*NSEG IN 2D, SO 2*NSEG IS ENOUGH TO STORE 3*NELEM
      flop1=>mesh%MSEG%X%R(        1:  nelem)
      flop2=>mesh%MSEG%X%R(  nelem+1:2*nelem)
      flop3=>mesh%MSEG%X%R(2*nelem+1:3*nelem)
      svol1=>     mesh%W%R(        1:  nelem)
      svol2=>     mesh%W%R(  nelem+1:2*nelem)
      svol3=>     mesh%W%R(2*nelem+1:3*nelem)
      dtlim1=>  mesh%M%X%R(        1:  nelem)
      dtlim2=>  mesh%M%X%R(  nelem+1:2*nelem)
      dtlim3=>  mesh%M%X%R(2*nelem+1:3*nelem)
!     ASSUMING THAT NSEG<3*NELEM
!     WILL BE DESTROYED WHEN DTLIM1 IS BUILT (BEWARE !!!!)
      fxmat=>mesh%M%X%R(1:nseg)
!
!-----------------------------------------------------------------------
!
!     VARIANTS (1 OR 2, DIFFERENT IMPLEMENTATION FOR SHARING VOLUMES
!               AT PREDICTOR LEVEL).
!
!     OPTPRE MUST BE EQUAL TO ITS VALUE IN POSITIVE_DEPTHS_ERIA !!!!!!!!
      optpre=1
!
!-----------------------------------------------------------------------
!
!     EXTRACTING OPTIONS, AND CONTROL
!
      iopt2=iopt/10
      IF(iopt2.NE.0.AND.iopt2.NE.1) THEN
        WRITE(lu,*) 'CVTRVF_ERIA: 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(h,t3)
      CALL cpstvc(h,t4)
      CALL cpstvc(h,t5)
      CALL cpstvc(h,t6)
      CALL cpstvc(h,t7)
!
!     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)
!     OPTIONAL COMPUTATION OF BOUNDARY FLUXES
      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
          limtra(i)=kddl
        ELSEIF(limtra(i).EQ.kddl.AND.t2%R(n).LT.0.d0) THEN
          limtra(i)=kdir
!         HERE FBOR MAY NOT HAVE BEEN GIVEN, FN TAKEN
          fbor%R(i)=fn%R(n)
        ELSE
!         INITIALISING FLBORTRA FOR POINTS NOT KDIR AND KDDL
          flbortra%R(i)=0.d0
        ENDIF
      ENDDO
!
!     INITIALIZING F AND INITIAL H
!
      CALL os('X=Y     ',x=f,y=fn)
      CALL os('X=Y     ',x=ht,y=hn)
!
!     FLUXES
!
      formul='HUGRADP         '
      IF(solsys.EQ.2) formul(8:8)='2'
!     COMPUTING THE 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
!     COPYING EBE FLUXES INTO FLOPOINT
      DO i=1,nelem
        flop1(i)=mesh%W%R(i)
        flop2(i)=mesh%W%R(i+nelmax)
        flop3(i)=mesh%W%R(i+2*nelmax)
      ENDDO
!     ASSEMBLING EBE FLUXES BY SEGMENT (TE1 SUIVI DE FALSE NON UTILISE)
!     HARDCODED OPTION 2 (MUST BE THE SAME AS IN POSITIVE_DEPTHS)
!     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.,2)
!     ASSEMBLING FXMAT IN PARALLEL
      IF(ncsize.GT.1) THEN
        CALL parcom2_seg(fxmat,fxmat,fxmat,mesh%NSEG,1,2,1,mesh,1,11)
      ENDIF
!     CHANGING FLUXES FROM POINTS INTO N FLUXES BETWEEN POINTS
      DO ielem = 1,nelem
        a1 = abs(flop1(ielem))
        a2 = abs(flop2(ielem))
        a3 = abs(flop3(ielem))
        IF(a1.GE.a2.AND.a1.GE.a3) THEN
!         ALL FLOW TO AND FROM NODE 1
          flop1(ielem)=-flop2(ielem)
          flop2(ielem)=0.d0
!         FLOP3(IELEM)= UNCHANGED!
        ELSEIF(a2.GE.a1.AND.a2.GE.a3) THEN
!         ALL FLOW TO AND FROM NODE 2
!         FLOP1(IELEM)= UNCHANGED!
          flop2(ielem)=-flop3(ielem)
          flop3(ielem)=0.d0
        ELSE
!         ALL FLOW TO AND FROM NODE 3
          flop3(ielem)=-flop1(ielem)
          flop1(ielem)=0.d0
!         FLOP2(IELEM)= UNCHANGED!
        ENDIF
      ENDDO
!
!     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 WILL BE AT LEAST C*DT
          ht%R(i)=ht%R(i)+c*dt
          IF(c.NE.0.d0) THEN
!                                           VALUE IN RAIN
            f%R(i)=f%R(i)+((c*dt)/ht%R(i))*(train-f%R(i))
          ENDIF
        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
!       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=mesh%NELEM
!
      DO i=1,remain
        indic_cpos(i)=i
      ENDDO
!
!     MAXIMUM INITIAL FLUX
!
      cprev=0.d0
      DO ir=1,nelem
        cprev=cprev+abs(flop1(ir))+abs(flop2(ir))+abs(flop3(ir))
      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
!
!     COMPUTING DEMAND (T7) AND OFFER (T1)
!
      CALL os('X=0     ',x=t7)
      CALL os('X=0     ',x=t1)
      IF(optpre.EQ.1) THEN
      DO ir=1,remain
        i=indic_cpos(ir)
        i1=mesh%IKLE%I(i        )
        i2=mesh%IKLE%I(i  +nelmax)
        i3=mesh%IKLE%I(i+2*nelmax)
!       A PRIORI AVAILABLE VOLUMES
        vol1=mesh%SURFAC%R(i)*ht%R(i1)*tiers
        vol2=mesh%SURFAC%R(i)*ht%R(i2)*tiers
        vol3=mesh%SURFAC%R(i)*ht%R(i3)*tiers
!       FLUXES LEAVING POINTS AGAIN
        f1= flop1(i)-flop3(i)
        f2=-flop1(i)+flop2(i)
        f3=-flop2(i)+flop3(i)
!       DEMAND OR OFFER
        IF(f1*dt.GT.vol1) THEN
!         POINT IN NEED OF EXTRA WATER
          t7%R(i1)=t7%R(i1)+f1*dt-vol1
        ELSE
!         POINT THAT CAN OFFER EXTRA WATER
!         BUT NOT MORE THAN ITS OWN VOLUME (AS WOULD BE THE CASE IF F1<0...)
          t1%R(i1)=t1%R(i1)+min(vol1,vol1-f1*dt)
        ENDIF
        IF(f2*dt.GT.vol2) THEN
          t7%R(i2)=t7%R(i2)+f2*dt-vol2
        ELSE
          t1%R(i2)=t1%R(i2)+min(vol2,vol2-f2*dt)
        ENDIF
        IF(f3*dt.GT.vol3) THEN
          t7%R(i3)=t7%R(i3)+f3*dt-vol3
        ELSE
          t1%R(i3)=t1%R(i3)+min(vol3,vol3-f3*dt)
        ENDIF
      ENDDO
      ELSEIF(optpre.EQ.2) THEN
      DO ir=1,remain
        i=indic_cpos(ir)
        i1=mesh%IKLE%I(i        )
        i2=mesh%IKLE%I(i  +nelmax)
        i3=mesh%IKLE%I(i+2*nelmax)
!       A PRIORI AVAILABLE VOLUMES
        vol1=mesh%SURFAC%R(i)*ht%R(i1)*tiers
        vol2=mesh%SURFAC%R(i)*ht%R(i2)*tiers
        vol3=mesh%SURFAC%R(i)*ht%R(i3)*tiers
!       FLUXES LEAVING POINTS AGAIN
        f1= flop1(i)-flop3(i)
        f2=-flop1(i)+flop2(i)
        f3=-flop2(i)+flop3(i)
!       DEMAND AND OFFER
        IF(f1.GT.0.d0) t7%R(i1)=t7%R(i1)+f1
        IF(f2.GT.0.d0) t7%R(i2)=t7%R(i2)+f2
        IF(f3.GT.0.d0) t7%R(i3)=t7%R(i3)+f3
        t1%R(i1)=t1%R(i1)+vol1
        t1%R(i2)=t1%R(i2)+vol2
        t1%R(i3)=t1%R(i3)+vol3
      ENDDO
      ENDIF
      IF(ncsize.GT.1) THEN
        CALL parcom(t7,2,mesh)
        CALL parcom(t1,2,mesh)
      ENDIF
!
!     PREPARING VOLUMES AND LIMITATION FOR PREDICTOR
!
      DO ir=1,remain
!
        i=indic_cpos(ir)
!
        i1=mesh%IKLE%I(i        )
        i2=mesh%IKLE%I(i  +nelmax)
        i3=mesh%IKLE%I(i+2*nelmax)
!
!       A PRIORI AVAILABLE VOLUMES FOR THIS ELEMENT
!
        f1= flop1(i)-flop3(i)
        f2=-flop1(i)+flop2(i)
        f3=-flop2(i)+flop3(i)
!
!       PRELIMINARY DISTRIBUTION OF VOLUMES BETWEEN ELEMENTS,
!       ACCORDING TO DEMAND AND OFFER
!
        IF(optpre.EQ.1) THEN
        vol1=mesh%SURFAC%R(i)*ht%R(i1)*tiers
        vol2=mesh%SURFAC%R(i)*ht%R(i2)*tiers
        vol3=mesh%SURFAC%R(i)*ht%R(i3)*tiers
        IF(f1*dt.GT.vol1) THEN
!         VOLUME TOO SMALL, TRYING TO INCREASE IT WITH GIFTS FROM OTHER
!         ELEMENTS
          IF(t7%R(i1).GT.t1%R(i1)) THEN
!           MORE DEMAND THAN OFFER, NEED TO SHARE
!                                  Commas very important to have T1/T7 <1
!                                  before multiplication
!                                  (                 )
            vol1=vol1+(f1*dt-vol1)*(t1%R(i1)/t7%R(i1))
          ELSE
!           MORE OFFER THAN DEMAND, THE POINT CAN HAVE THE REQUESTED VOLUME
            vol1=f1*dt
          ENDIF
        ELSE
!         VOLUME LARGE ENOUGH, GIVING WATER TO OTHER ELEMENTS ?
          IF(t1%R(i1).GT.t7%R(i1)) THEN
!           MORE OFFER THAN DEMAND, GIVING ENOUGH, BUT NOT ALL
            vol1=vol1-min(vol1,vol1-f1*dt)*(t7%R(i1)/t1%R(i1))
          ELSE
!           MORE DEMAND THAN OFFER, KEEPING ONLY WHAT IS NECESSARY
            vol1=max(f1,0.d0)*dt
          ENDIF
        ENDIF
        IF(f2*dt.GT.vol2) THEN
          IF(t7%R(i2).GT.t1%R(i2)) THEN
            vol2=vol2+(f2*dt-vol2)*(t1%R(i2)/t7%R(i2))
          ELSE
            vol2=f2*dt
          ENDIF
        ELSE
          IF(t1%R(i2).GT.t7%R(i2)) THEN
            vol2=vol2-min(vol2,vol2-f2*dt)*(t7%R(i2)/t1%R(i2))
          ELSE
            vol2=max(f2,0.d0)*dt
          ENDIF
        ENDIF
        IF(f3*dt.GT.vol3) THEN
          IF(t7%R(i3).GT.t1%R(i3)) THEN
            vol3=vol3+(f3*dt-vol3)*(t1%R(i3)/t7%R(i3))
          ELSE
            vol3=f3*dt
          ENDIF
        ELSE
          IF(t1%R(i3).GT.t7%R(i3)) THEN
            vol3=vol3-min(vol3,vol3-f3*dt)*(t7%R(i3)/t1%R(i3))
          ELSE
            vol3=max(f3,0.d0)*dt
          ENDIF
        ENDIF
        ELSEIF(optpre.EQ.2) THEN
        IF(t7%R(i1).GT.1.d-30) THEN
!                       ( THIS IS IMPORTANT
          vol1=t1%R(i1)*(max(f1,0.d0)/t7%R(i1))
        ELSE
          vol1=mesh%SURFAC%R(i)*h%R(i1)*tiers
        ENDIF
        IF(t7%R(i2).GT.1.d-30) THEN
          vol2=t1%R(i2)*(max(f2,0.d0)/t7%R(i2))
        ELSE
          vol2=mesh%SURFAC%R(i)*h%R(i2)*tiers
        ENDIF
        IF(t7%R(i3).GT.1.d-30) THEN
          vol3=t1%R(i3)*(max(f3,0.d0)/t7%R(i3))
        ELSE
          vol3=mesh%SURFAC%R(i)*h%R(i3)*tiers
        ENDIF
        ENDIF
!
!       SAVING VOLUMES FOR CORRECTOR (ACTUALLY NOT USED BY PREDICTOR, THOUGH
!                                     DONE FOR IT, BUT IT GIVES THE LIMITATIONS...)
!
        svol1(i)=vol1
        svol2(i)=vol2
        svol3(i)=vol3
!
!       LIMITATION OF FLUXES, FIRST IN TERMS OF LIMITED TIME-STEP
!
          IF(f1*dt.GT.vol1) THEN
            dt1=dt*(vol1/(f1*dt))
          ELSE
            dt1=dt
          ENDIF
          IF(f2*dt.GT.vol2) THEN
            dt2=dt*(vol2/(f2*dt))
          ELSE
            dt2=dt
          ENDIF
          IF(f3*dt.GT.vol3) THEN
            dt3=dt*(vol3/(f3*dt))
          ELSE
            dt3=dt
          ENDIF
!
!       NOW LIMITED VOLUMES TRANSITING BETWEEN POINTS (1/DT MISSING)
!       A SEGMENT GETS THE STRONGEST LIMITATION OF ITS TWO POINTS
!       THIS WILL CREATE NO NEGATIVE VOLUME, AS IT COULD BE FEARED,
!       BECAUSE WE HAVE ONLY THE 1-TARGET CASE AND THE 2-TARGET CASE
!       THUS ALL FLUXES LEAVING OR ARRIVING TO A POINT HAVE THE SAME
!       SIGN OR ARE ZERO !
!
        dtlim1(i)=min(dt1,dt2)
        dtlim2(i)=min(dt2,dt3)
        dtlim3(i)=min(dt3,dt1)
!
      ENDDO
!
!
!     NOW SEQUENCE PREDICTOR-CORRECTOR
!
!
!     INITIALISING T3, FUTURE VALUE OF F BECAUSE POINTS THAT DO NOT RECEIVE
!                      WATER WILL NOT BE TREATED
!     T4 IS THE EVOLUTION OF VOLUME OF WATER, HERE INITIALISED TO 0
!     T5 IS THE EVOLUTION OF MASSES OF TRACER, HERE INITIALISED TO 0
!     T6 IS THE INITIAL TRACER*DEPTH
      CALL os('X=0     ',x=t4)
      CALL os('X=0     ',x=t5)
      DO i=1,f%DIM1
        t6%R(i)=f%R(i)*ht%R(i)
      ENDDO
      DO i=1,f%DIM1
        t3%R(i)=f%R(i)
      ENDDO
!
!     COMPUTING MIN AND MAX FOR FINAL CLIPPING, THIS WILL ENSURE A GLOBAL
!     MONOTONICITY, THAT COULD OTHERWISE BE VIOLATED BY TRUNCATION ERRORS
!     THIS COULD BE DONE LOCALLY, BUT IT WOULD BE MORE EXPENSIVE
!     FOR LOCAL MIN AND MAX SEE CVTRVF.F
!
      fmin=f%R(1)
      fmax=f%R(1)
      DO i=2,f%DIM1
        fmin=min(fmin,f%R(i))
        fmax=max(fmax,f%R(i))
      ENDDO
      IF(ncsize.GT.1) THEN
        fmin=p_min(fmin)
        fmax=p_max(fmax)
      ENDIF
!
!     PREDICTOR
!
      DO ir=1,remain
!
        i=indic_cpos(ir)
!
        i1=mesh%IKLE%I(i        )
        i2=mesh%IKLE%I(i  +nelmax)
        i3=mesh%IKLE%I(i+2*nelmax)
!
        fp1=flop1(i)*dtlim1(i)
        fp2=flop2(i)*dtlim2(i)
        fp3=flop3(i)*dtlim3(i)
!
!       CORRESPONDING VARIATIONS OF VOLUMES OF POINTS (DT MISSING SO OK)
!
        t4%R(i1)=t4%R(i1)-( fp1-fp3)
        t4%R(i2)=t4%R(i2)-(-fp1+fp2)
        t4%R(i3)=t4%R(i3)-(-fp2+fp3)
!
!       VARIATIONS OF MASSES ADDED AND REMOVED IN REDUCED FORM
!
        t5%R(i1)=t5%R(i1)-( fp1-fp3)*f%R(i1)
        t5%R(i2)=t5%R(i2)-(-fp1+fp2)*f%R(i2)
        t5%R(i3)=t5%R(i3)-(-fp2+fp3)*f%R(i3)
!
        fi1=-min(fp1,0.d0)*(f%R(i2)-f%R(i1))
     &      +max(fp3,0.d0)*(f%R(i3)-f%R(i1))
        fi2=+max(fp1,0.d0)*(f%R(i1)-f%R(i2))
     &      -min(fp2,0.d0)*(f%R(i3)-f%R(i2))
        fi3=+max(fp2,0.d0)*(f%R(i2)-f%R(i3))
     &      -min(fp3,0.d0)*(f%R(i1)-f%R(i3))
!
!       PSI LIMITATION
!
        fitot=fi1+fi2+fi3
!
        IF(fitot.GT.eps_flux) THEN
!         PSI REDUCTION
          beta1=max(fi1,0.d0)
          beta2=max(fi2,0.d0)
          beta3=max(fi3,0.d0)
          coef=fitot/(beta1+beta2+beta3)
          fi1=beta1*coef
          fi2=beta2*coef
          fi3=beta3*coef
        ELSEIF(fitot.LT.-eps_flux) THEN
!         PSI REDUCTION
          beta1=min(fi1,0.d0)
          beta2=min(fi2,0.d0)
          beta3=min(fi3,0.d0)
          coef=fitot/(beta1+beta2+beta3)
          fi1=beta1*coef
          fi2=beta2*coef
          fi3=beta3*coef
        ELSE
!         NO REDUCTION
        ENDIF
!
        t5%R(i1)=t5%R(i1)+fi1
        t5%R(i2)=t5%R(i2)+fi2
        t5%R(i3)=t5%R(i3)+fi3
!
      ENDDO
!
!     ADDING THE EVOLUTIONS TO THE DEPTHS AND TRACERS
!     AFTER ASSEMBLY AT INTERFACES AND AFTER
!     CHANGING VOLUMES INTO DEPTHS.
!
      IF(ncsize.GT.1) THEN
        CALL parcom(t4,2,mesh)
        CALL parcom(t5,2,mesh)
      ENDIF
      CALL os('X=XY    ',x=t4,y=unsv2d)
      CALL os('X=XY    ',x=t5,y=unsv2d)
!     UPDATING DEPOTH AND TRACER
      DO i=1,f%DIM1
!       ADDING INITIAL QUANTITY T6 AND VARIATION T5
        t6%R(i)=t6%R(i)+t5%R(i)
!       NOW T5 USED TO STORE OLD DEPTH
        t5%R(i)=ht%R(i)
!       NEW DEPTH
        ht%R(i)=ht%R(i)+t4%R(i)
        IF(ht%R(i).GT.1.d-15) THEN
!         NEW VALUE OF TRACER
          t3%R(i)=t6%R(i)/ht%R(i)
        ENDIF
!       TO COPE WITH TRUNCATION ERRORS AND AVOID NEGATIVE VALUES
!       IF HT NEGATIVE FOR OTHER REASONS IT WILL MAKE MASS ERRORS
        ht%R(i)=max(0.d0,ht%R(i))
!       ENSURING A GLOBAL MONOTONICITY OF T3
        t3%R(i)=min(fmax,max(fmin,t3%R(i)))
      ENDDO
!
!     CORRECTOR
!
      IF(nco_dist.GT.0) THEN
!
        DO i=1,f%DIM1
          t7%R(i)=0.d0
          t1%R(i)=0.d0
        ENDDO
!
!       EVALUATING OFFER (T1) AND DEMAND (T7)
!
        DO ir=1,remain
          i=indic_cpos(ir)
!
          i1=mesh%IKLE%I(i        )
          i2=mesh%IKLE%I(i  +nelmax)
          i3=mesh%IKLE%I(i+2*nelmax)
!         LIMITED VOLUMES TRAVELLING BETWEEN POINTS
          fp1=flop1(i)*dtlim1(i)
          fp2=flop2(i)*dtlim2(i)
          fp3=flop3(i)*dtlim3(i)
!         VOLUMES OF BASIC SOLUTION AT T(N+1)
          vol1=svol1(i)-( fp1-fp3)
          vol2=svol2(i)-(-fp1+fp2)
          vol3=svol3(i)-(-fp2+fp3)
          IF(vol1.LT.mesh%SURFAC%R(i)*ht%R(i1)*tiers) THEN
            t7%R(i1)=t7%R(i1)+mesh%SURFAC%R(i)*ht%R(i1)*tiers-vol1
          ELSE
            t1%R(i1)=t1%R(i1)+vol1+min( fp1,0.d0)+min(-fp3,0.d0)
          ENDIF
          IF(vol2.LT.mesh%SURFAC%R(i)*ht%R(i2)*tiers) THEN
            t7%R(i2)=t7%R(i2)+mesh%SURFAC%R(i)*ht%R(i2)*tiers-vol2
          ELSE
            t1%R(i2)=t1%R(i2)+vol2+min(-fp1,0.d0)+min( fp2,0.d0)
          ENDIF
          IF(vol3.LT.mesh%SURFAC%R(i)*ht%R(i3)*tiers) THEN
            t7%R(i3)=t7%R(i3)+mesh%SURFAC%R(i)*ht%R(i3)*tiers-vol3
          ELSE
            t1%R(i3)=t1%R(i3)+vol3+min(-fp2,0.d0)+min( fp3,0.d0)
          ENDIF
        ENDDO
!
        IF(ncsize.GT.1) THEN
!         GLOBAL OFFER (T1) AND DEMAND (T7)
          CALL parcom(t1,2,mesh)
          CALL parcom(t7,2,mesh)
        ENDIF
!
!       VOLUMES THAT WILL BE USED FOR THE DERIVATIVE IN TIME
!
        DO ir=1,remain
!
          i=indic_cpos(ir)
!
          i1=mesh%IKLE%I(i        )
          i2=mesh%IKLE%I(i  +nelmax)
          i3=mesh%IKLE%I(i+2*nelmax)
!
!         LIMITED VOLUMES BETWEEN POINTS
!
          fp1=flop1(i)*dtlim1(i)
          fp2=flop2(i)*dtlim2(i)
          fp3=flop3(i)*dtlim3(i)
!         LOCAL VOLUMES AFTER PREDICTOR
          vol1=svol1(i)-( fp1-fp3)
          vol2=svol2(i)-(-fp1+fp2)
          vol3=svol3(i)-(-fp2+fp3)
!         VOLUME1
          IF(vol1.LT.mesh%SURFAC%R(i)*ht%R(i1)*tiers) THEN
            IF(t7%R(i1).GT.t1%R(i1)) THEN
!             MORE DEMAND THAN OFFER, VOL1 GETS ONLY A SHARE
              vol1=vol1+(mesh%SURFAC%R(i)*ht%R(i1)*tiers-vol1)
     &                 *(t1%R(i1)/t7%R(i1))
            ELSE
!             VOL1 GETS ALL
              vol1=mesh%SURFAC%R(i)*ht%R(i1)*tiers
            ENDIF
          ELSE
            IF(t1%R(i1).GT.t7%R(i1)) THEN
!             MORE OFFER THAN DEMAND, VOL1 GIVES A SHARE ONLY
              vol1=vol1-(vol1+(min( fp1,0.d0)+min(-fp3,0.d0)))
     &                 *(t7%R(i1)/t1%R(i1))
            ELSE
!             MORE DEMAND THAN OFFER, VOL1 GIVES ALL IT CAN
              vol1=-(min( fp1,0.d0)+min(-fp3,0.d0))
            ENDIF
          ENDIF
!         VOLUME 2
          IF(vol2.LT.mesh%SURFAC%R(i)*ht%R(i2)*tiers) THEN
            IF(t7%R(i2).GT.t1%R(i2)) THEN
!             MORE DEMAND THAN OFFER, VOL2 GETS ONLY A SHARE
              vol2=vol2+(mesh%SURFAC%R(i)*ht%R(i2)*tiers-vol2)
     &                 *(t1%R(i2)/t7%R(i2))
            ELSE
!             VOL2 GETS ALL
              vol2=mesh%SURFAC%R(i)*ht%R(i2)*tiers
            ENDIF
          ELSE
            IF(t1%R(i2).GT.t7%R(i2)) THEN
!             MORE OFFER THAN DEMAND, VOL2 GIVES A SHARE ONLY
              vol2=vol2-(vol2+(min(-fp1,0.d0)+min( fp2,0.d0)))
     &                 *(t7%R(i2)/t1%R(i2))
            ELSE
!             MORE DEMAND THAN OFFER, VOL1 GIVES ALL IT CAN
              vol2=-(min(-fp1,0.d0)+min( fp2,0.d0))
            ENDIF
          ENDIF
!         VOLUME 3
          IF(vol3.LT.mesh%SURFAC%R(i)*ht%R(i3)*tiers) THEN
            IF(t7%R(i3).GT.t1%R(i3)) THEN
!             MORE DEMAND THAN OFFER, VOL1 GETS ONLY A SHARE
              vol3=vol3+(mesh%SURFAC%R(i)*ht%R(i3)*tiers-vol3)
     &                 *(t1%R(i3)/t7%R(i3))
            ELSE
!             VOL3 GETS ALL
              vol3=mesh%SURFAC%R(i)*ht%R(i3)*tiers
            ENDIF
          ELSE
            IF(t1%R(i3).GT.t7%R(i3)) THEN
!             MORE OFFER THAN DEMAND, VOL3 GIVES A SHARE ONLY
              vol3=vol3-(vol3+(min(-fp2,0.d0)+min( fp3,0.d0)))
     &                 *(t7%R(i3)/t1%R(i3))
            ELSE
!             MORE DEMAND THAN OFFER, VOL3 GIVES ALL IT CAN
              vol3=-(min(-fp2,0.d0)+min( fp3,0.d0))
            ENDIF
          ENDIF
          svol1(i)=vol1
          svol2(i)=vol2
          svol3(i)=vol3
        ENDDO
!
        DO n=1,nco_dist
!
          DO i=1,f%DIM1
!           INITIAL QUANTITY
            t6%R(i)=f%R(i)*ht%R(i)
          ENDDO
!
          IF(optadv.EQ.2) THEN
!           FOR SECOND ORDER, COMPUTING POSSIBLE MIN AND MAX
!           MIN IN T1 AND MAX IN T7
            DO i=1,f%DIM1
              t1%R(i)=f%R(i)
              t7%R(i)=f%R(i)
            ENDDO
            DO ir=1,remain
              i=indic_cpos(ir)
              i1=mesh%IKLE%I(i        )
              i2=mesh%IKLE%I(i  +nelmax)
              i3=mesh%IKLE%I(i+2*nelmax)
              localmin=min( f%R(i1), f%R(i2), f%R(i3),
     &                     t3%R(i1),t3%R(i2),t3%R(i3))
              localmax=max( f%R(i1), f%R(i2), f%R(i3),
     &                     t3%R(i1),t3%R(i2),t3%R(i3))
              t1%R(i1)=min(t1%R(i1),localmin)
              t1%R(i2)=min(t1%R(i2),localmin)
              t1%R(i3)=min(t1%R(i3),localmin)
              t7%R(i1)=max(t7%R(i1),localmax)
              t7%R(i2)=max(t7%R(i2),localmax)
              t7%R(i3)=max(t7%R(i3),localmax)
            ENDDO
            IF(ncsize.GT.1) THEN
!                            4:MIN
              CALL parcom(t1,4,mesh)
!                            3:MAX
              CALL parcom(t7,3,mesh)
            ENDIF
          ENDIF
!
          CALL os('X=0     ',x=t5)
!
!         NOW THE REAL CORRECTOR
!
          IF(optadv.EQ.1) THEN
!
!           FIRST ORDER IN TIME
!
            DO ir=1,remain
!
              i=indic_cpos(ir)
              i1=mesh%IKLE%I(i        )
              i2=mesh%IKLE%I(i  +nelmax)
              i3=mesh%IKLE%I(i+2*nelmax)
!
!             LIMITED VOLUMES BETWEEN POINTS
!
              fp1=flop1(i)*dtlim1(i)
              fp2=flop2(i)*dtlim2(i)
              fp3=flop3(i)*dtlim3(i)
!
              fi1=(t3%R(i1)-f%R(i1))*svol1(i)
              fi2=(t3%R(i2)-f%R(i2))*svol2(i)
              fi3=(t3%R(i3)-f%R(i3))*svol3(i)
!
!             ADDING THE DERIVATIVE THAT WILL BE REMOVED IN UPWIND FORM
!             THIS IS DONE ONLY FOR REMAINING ELEMENTS.
!
              t5%R(i1)=t5%R(i1)+fi1
              t5%R(i2)=t5%R(i2)+fi2
              t5%R(i3)=t5%R(i3)+fi3
!
!             THEN FLUXES
!
              fip1=-min(fp1,0.d0)*(f%R(i2)-f%R(i1))
     &             +max(fp3,0.d0)*(f%R(i3)-f%R(i1))
              fip2=+max(fp1,0.d0)*(f%R(i1)-f%R(i2))
     &             -min(fp2,0.d0)*(f%R(i3)-f%R(i2))
              fip3=+max(fp2,0.d0)*(f%R(i2)-f%R(i3))
     &             -min(fp3,0.d0)*(f%R(i1)-f%R(i3))
!
!             THIS EXTRA PSI REDUCTION IS NOT NECESSARY
!             BUT GIVES SLIGHTLY BETTER RESULTS
              fitot=fip1+fip2+fip3
              IF(fitot.GT.eps_flux) THEN
!               PSI REDUCTION
                beta1=max(fip1,0.d0)
                beta2=max(fip2,0.d0)
                beta3=max(fip3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fip1=beta1*coef
                fip2=beta2*coef
                fip3=beta3*coef
              ELSEIF(fitot.LT.-eps_flux) THEN
!               PSI REDUCTION
                beta1=min(fip1,0.d0)
                beta2=min(fip2,0.d0)
                beta3=min(fip3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fip1=beta1*coef
                fip2=beta2*coef
                fip3=beta3*coef
              ELSE
!               NO REDUCTION
              ENDIF
!
!             ADDING TO FINAL CONTRIBUTIONS THAT WILL BE REDUCED AGAIN
!
              fi1=-fi1+fip1
              fi2=-fi2+fip2
              fi3=-fi3+fip3
!
!             PSI LIMITATION
!
              fitot=fi1+fi2+fi3
!
              IF(fitot.GT.eps_flux) THEN
!               PSI REDUCTION
                beta1=max(fi1,0.d0)
                beta2=max(fi2,0.d0)
                beta3=max(fi3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fi1=beta1*coef
                fi2=beta2*coef
                fi3=beta3*coef
              ELSEIF(fitot.LT.-eps_flux) THEN
!               PSI REDUCTION
                beta1=min(fi1,0.d0)
                beta2=min(fi2,0.d0)
                beta3=min(fi3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fi1=beta1*coef
                fi2=beta2*coef
                fi3=beta3*coef
              ELSE
!               NO REDUCTION
              ENDIF
!
              t5%R(i1)=t5%R(i1)+fi1
              t5%R(i2)=t5%R(i2)+fi2
              t5%R(i3)=t5%R(i3)+fi3
!
            ENDDO
!
          ELSEIF(optadv.EQ.2) THEN
!
!           SECOND ORDER IN TIME
!
            DO ir=1,remain
!
              i=indic_cpos(ir)
              i1=mesh%IKLE%I(i        )
              i2=mesh%IKLE%I(i  +nelmax)
              i3=mesh%IKLE%I(i+2*nelmax)
!
!             LIMITED VOLUMES BETWEEN POINTS
!
              fp1=flop1(i)*dtlim1(i)
              fp2=flop2(i)*dtlim2(i)
              fp3=flop3(i)*dtlim3(i)
!
              fi1=(f%R(i1)-t3%R(i1))*svol1(i)
              fi2=(f%R(i2)-t3%R(i2))*svol2(i)
              fi3=(f%R(i3)-t3%R(i3))*svol3(i)
!
!             ADDING THE DERIVATIVE THAT WILL BE REMOVED IN UPWIND FORM
!
              t5%R(i1)=t5%R(i1)-fi1
              t5%R(i2)=t5%R(i2)-fi2
              t5%R(i3)=t5%R(i3)-fi3
!
!             THEN FLUXES
!
              vold1=max(0.d0,svol1(i)+min(fp1,0.d0)-max(fp3,0.d0))
              vold2=max(0.d0,svol2(i)-max(fp1,0.d0)+min(fp2,0.d0))
              vold3=max(0.d0,svol3(i)-max(fp2,0.d0)+min(fp3,0.d0))
              a1=2.d0*vold1/max(1.d-30, max(fp1,0.d0)-min(fp3,0.d0))
              a2=2.d0*vold2/max(1.d-30,-min(fp1,0.d0)+max(fp2,0.d0))
              a3=2.d0*vold3/max(1.d-30,-min(fp2,0.d0)+max(fp3,0.d0))
              fs1=min(t3%R(i1),f%R(i1)+a1*(f%R(i1)-t1%R(i1)))
              fs1=max(fs1     ,f%R(i1)+a1*(f%R(i1)-t7%R(i1)))
              fs2=min(t3%R(i2),f%R(i2)+a2*(f%R(i2)-t1%R(i2)))
              fs2=max(fs2     ,f%R(i2)+a2*(f%R(i2)-t7%R(i2)))
              fs3=min(t3%R(i3),f%R(i3)+a3*(f%R(i3)-t1%R(i3)))
              fs3=max(fs3     ,f%R(i3)+a3*(f%R(i3)-t7%R(i3)))
!
              fip1=0.5d0*(-min(fp1,0.d0)*(fs2-f%R(i1)+f%R(i2)-f%R(i1))
     &                    +max(fp3,0.d0)*(fs3-f%R(i1)+f%R(i3)-f%R(i1))
     &                    -max(fp1,0.d0)*(fs1-f%R(i1))
     &                    +min(fp3,0.d0)*(fs1-f%R(i1)))
              fip2=0.5d0*(+max(fp1,0.d0)*(fs1-f%R(i2)+f%R(i1)-f%R(i2))
     &                    -min(fp2,0.d0)*(fs3-f%R(i2)+f%R(i3)-f%R(i2))
     &                    +min(fp1,0.d0)*(fs2-f%R(i2))
     &                    -max(fp2,0.d0)*(fs2-f%R(i2)))
              fip3=0.5d0*(+max(fp2,0.d0)*(fs2-f%R(i3)+f%R(i2)-f%R(i3))
     &                    -min(fp3,0.d0)*(fs1-f%R(i3)+f%R(i1)-f%R(i3))
     &                    +min(fp2,0.d0)*(fs3-f%R(i3))
     &                    -max(fp3,0.d0)*(fs3-f%R(i3)))
!
!             GO TO 2000
!             THIS EXTRA PSI REDUCTION IS NOT NECESSARY
!             BUT GIVES SLIGHTLY BETTER RESULTS
              fitot=fip1+fip2+fip3
              IF(fitot.GT.eps_flux) THEN
!               PSI REDUCTION
                beta1=max(fip1,0.d0)
                beta2=max(fip2,0.d0)
                beta3=max(fip3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fip1=beta1*coef
                fip2=beta2*coef
                fip3=beta3*coef
              ELSEIF(fitot.LT.-eps_flux) THEN
!               PSI REDUCTION
                beta1=min(fip1,0.d0)
                beta2=min(fip2,0.d0)
                beta3=min(fip3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fip1=beta1*coef
                fip2=beta2*coef
                fip3=beta3*coef
              ELSE
!               NO REDUCTION
              ENDIF
!2000          CONTINUE
!
!             ADDING TO FINAL CONTRIBUTIONS THAT WILL BE REDUCED AGAIN
!
              fi1=fi1+fip1
              fi2=fi2+fip2
              fi3=fi3+fip3
!
!             PSI LIMITATION
!
              fitot=fi1+fi2+fi3
!
              IF(fitot.GT.eps_flux) THEN
!               PSI REDUCTION
                beta1=max(fi1,0.d0)
                beta2=max(fi2,0.d0)
                beta3=max(fi3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fi1=beta1*coef
                fi2=beta2*coef
                fi3=beta3*coef
              ELSEIF(fitot.LT.-eps_flux) THEN
!               PSI REDUCTION
                beta1=min(fi1,0.d0)
                beta2=min(fi2,0.d0)
                beta3=min(fi3,0.d0)
                coef=fitot/(beta1+beta2+beta3)
                fi1=beta1*coef
                fi2=beta2*coef
                fi3=beta3*coef
              ELSE
!               NO REDUCTION
              ENDIF
!
              t5%R(i1)=t5%R(i1)+fi1
              t5%R(i2)=t5%R(i2)+fi2
              t5%R(i3)=t5%R(i3)+fi3
!
            ENDDO
!
          ELSE
            WRITE(lu,*) 'UNKNOWN OPTION IN CVTRVF_ERIA: ',optadv
            CALL plante(1)
            stop
          ENDIF
!
!         ADDING THE EVOLUTIONS TO THE TRACERS
!         AFTER ASSEMBLY AT INTERFACES AND AFTER
!         CHANGING VOLUMES INTO DEPTHS.
!
          IF(ncsize.GT.1) CALL parcom(t5,2,mesh)
          CALL os('X=XY    ',x=t5,y=unsv2d)
!         NEW AVERAGE VALUE OF TRACER
          DO i=1,f%DIM1
            IF(ht%R(i).GT.1.d-15) THEN
              t3%R(i)=(t6%R(i)+t5%R(i))/ht%R(i)
!             ENSURING A GLOBAL MONOTONICITY
              t3%R(i)=min(fmax,max(fmin,t3%R(i)))
            ENDIF
          ENDDO
!
!         END OF THE LOOP OF CORRECTIONS
!
        ENDDO
      ENDIF
!
!     SETTING THE FINAL VALUE OF THE TRACER FOR THIS SUB-ITERATION
      DO i=1,f%DIM1
        f%R(i)=t3%R(i)
      ENDDO
!
!     IF REMAINING FLUXES, THE ELEMENT IS KEPT IN THE LIST
      newremain=0
      c=0.d0
!
      DO ir=1,remain
        i=indic_cpos(ir)
        IF(dtlim1(i).EQ.dt.AND.dtlim2(i).EQ.dt.AND.
     &     dtlim3(i).EQ.dt) THEN
          flop1(i)=0.d0
          flop2(i)=0.d0
          flop3(i)=0.d0
        ELSE
          newremain=newremain+1
!         BEFORE NEWREMAIN: FOR NEXT ITERATION
!         AFTER  NEWREMAIN: STILL VALID FOR NEXT ITERATION
          indic_cpos(newremain)=i
          flop1(i)=flop1(i)*(1.d0-dtlim1(i)*surdt)
          flop2(i)=flop2(i)*(1.d0-dtlim2(i)*surdt)
          flop3(i)=flop3(i)*(1.d0-dtlim3(i)*surdt)
          c=c+abs(flop1(i))+abs(flop2(i))+abs(flop3(i))
        ENDIF
      ENDDO
!
      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)+c*dt
!         MONOTONICITY NON OBEYED HERE OF COURSE
!                        VALUE IN VAPOR
!         F%R(I)=F%R(I)+(0.D0-F%R(I))*((C*DT)/MAX(HT%R(I),1.D-4))
          f%R(i)=f%R(i)*(1.d0-((c*dt)/max(ht%R(i),1.d-4)))
        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
!
      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)-DT*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,102) niter
      ENDIF
!
      IF(niter.EQ.nitmax) THEN
        WRITE(lu,103) niter
      ENDIF
!
102   FORMAT(' CVTRVF_ERIA (SCHEME ERIA, 15): ',1i3,' ITERATIONS')
103   FORMAT(' CVTRVF_ERIA (SCHEME ERIA, 15): ',1i3,' ITERATIONS',
     &  ' = MAXIMUM',
     &  /,'INCREASE MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES')
!
!-----------------------------------------------------------------------
!
      RETURN
      END