flux_ef_vf_3.f

Go to the documentation of this file.

!                   ***********************
                    SUBROUTINE flux_ef_vf_3
!                   ***********************
!
     &(phiel,nelem,nelmax,eltseg,oriseg,fxmatpar,nseg,ikle,npoin,
     & fn,fi_i,su,hdfdt,teta,yaflulim,flulim,yaflulimebe,flulimebe)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    Equivalent of FLUX_EF_VF, but only for PSI scheme, and the
!+        result is given in terms of contribution per point, not fluxes
!+        between points, and takes a derivative in time into account.
!
!warning  When both YAFLULIM and YAFLULIMEBE are true only the latter is
!+        taken into account !
!
!history  J-M HERVOUET & SARA PAVAN (EDF LAB, LNHE)
!+        04/08/2015
!+        V7P1
!+   First version. The N and PSI versions are not different here.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        05/11/2016
!+        V7P3
!+   Now two options of flux limitation offered.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        09/09/2017
!+        V7P3
!+   Argument NELMAX added for correctly dimensioning arrays.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ELTSEG         |-->| GIVES THE SEGMENTS IN A TRIANGLE
!| FI_I           |<--| CONTRIBUTIONS TO POINTS
!| FLULIM         |-->| LIMITATION OF FLUXES BY SEGMENT
!| FLULIMEBE      |-->| LIMITATION OF FLUXES BY SEGMENT BUT EBE
!| FN             |-->| TRACER AT TIME T(N)
!| FXMARPAR       |-->| FLUXES ASSEMBLED IN //
!| HDFDT          |-->| DEPTH*DH/DT
!| IKLE           |-->| CONNECTIVITY TABLE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELEM          |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| NSEG           |-->| NUMBER OF SEGMENTS
!| ORISEG         |-->| GIVES THE ORIENTATION OF SEGMENTS IN A TRIANGLE
!| PHIEL          |-->| PER ELEMENT, FLUXES LEAVING POINTS
!| SU             |-->| SURFACE OF TRIANGLES
!| TETA           |-->| LOCAL IMPLICITATION
!| YAFLULIM       |-->| IF YES, A LIMITATION OF FLUXES IS GIVEN BY SEGMENT
!| YAFLULIMEBE    |-->| IF YES, AN EBE LIMITATION OF FLUXES IS GIVEN
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      USE bief, ex_flux_ef_vf_3 => flux_ef_vf_3
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NELEM,NELMAX,NPOIN,NSEG
      INTEGER, INTENT(IN)             :: IKLE(nelmax,3)
      INTEGER, INTENT(IN)             :: ELTSEG(nelmax,3)
      INTEGER, INTENT(IN)             :: ORISEG(nelmax,3)
      LOGICAL, INTENT(IN)             :: YAFLULIM,YAFLULIMEBE
      DOUBLE PRECISION, INTENT(IN)    :: PHIEL(nelmax,3),TETA(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: FLULIMEBE(*)
      DOUBLE PRECISION, INTENT(IN)    :: FXMATPAR(nseg),FLULIM(*)
      DOUBLE PRECISION, INTENT(INOUT) :: FI_I(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: SU(nelem),HDFDT(npoin)
      TYPE(bief_obj), INTENT(IN)      :: FN
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,I1,I2,I3,I,ISEG1,ISEG2,ISEG3
      DOUBLE PRECISION K1,K2,K3,FN1,FN2,FN3,BETA1,BETA2,BETA3
      DOUBLE PRECISION FINCORR1,FINCORR2,FINCORR3,PHITCOR,COEF
!     PSI FLUXES WITH FN
      DOUBLE PRECISION FP21,FP32,FP13,FP12,FP23,FP31
      DOUBLE PRECISION MIN12,MIN23,MIN13,MT1,MT2,MT3
!
      INTRINSIC min,max
!
      DOUBLE PRECISION, PARAMETER :: EPSPHI=1.d-12
      DOUBLE PRECISION, PARAMETER :: TIERS =1.d0/3.d0
!
!-----------------------------------------------------------------------
!
      DO i=1,npoin
        fi_i(i)=0.d0
      ENDDO
!
!     DEPENDING ON THE FLUX LIMITATION...
!
      IF(yaflulimebe) THEN
!
!       AS WHEN YAFLULIM=.TRUE., BUT WITH FLULIM ASSUMED TO BE 1
!
        DO ielem = 1,nelem
!
          k1 = -phiel(ielem,1)
          k2 = -phiel(ielem,2)
          k3 = -phiel(ielem,3)
!
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
!
          iseg1=eltseg(ielem,1)
          iseg2=eltseg(ielem,2)
          iseg3=eltseg(ielem,3)
!
!         STARTS WITH N-SCHEME (EQUIVALENT TO LEO POSTMA'S IMPLEMENTATION)
!         FIJ HERE LIKE LAMBDA(I,J) IN BOOK
!         I.E. FIJ>0 MEANS FLUX FROM J TO I !!!!!
!
          fp12=max(min(k1,-k2),0.d0)*flulimebe((1-1)*nelmax+ielem)
          fp23=max(min(k2,-k3),0.d0)*flulimebe((2-1)*nelmax+ielem)
          fp31=max(min(k3,-k1),0.d0)*flulimebe((3-1)*nelmax+ielem)
          fp21=max(min(k2,-k1),0.d0)*flulimebe((1-1)*nelmax+ielem)
          fp32=max(min(k3,-k2),0.d0)*flulimebe((2-1)*nelmax+ielem)
          fp13=max(min(k1,-k3),0.d0)*flulimebe((3-1)*nelmax+ielem)
!
!         CORRECTING THE FLUXES WHEN THEIR SIGN IS NOT THE SAME
!         AS THE ASSEMBLED VALUE, KNOWING THAT ALL THE FPIJ ARE
!         POSITIVE BY CONSTRUCTION
!
!         SEGMENT 1
!
          IF(oriseg(ielem,1).EQ.1) THEN
            IF(fxmatpar(iseg1).GT.0.d0) THEN
              fp12=0.d0
              IF(fp21.GT. fxmatpar(iseg1)) fp21=fxmatpar(iseg1)
            ELSE
              fp21=0.d0
              IF(fp12.GT.-fxmatpar(iseg1)) fp12=-fxmatpar(iseg1)
            ENDIF
          ELSE
            IF(fxmatpar(iseg1).GT.0.d0) THEN
              fp21=0.d0
              IF(fp12.GT. fxmatpar(iseg1)) fp12=fxmatpar(iseg1)
            ELSE
              fp12=0.d0
              IF(fp21.GT.-fxmatpar(iseg1)) fp21=-fxmatpar(iseg1)
            ENDIF
          ENDIF
!
!         SEGMENT 2
!
          IF(oriseg(ielem,2).EQ.1) THEN
            IF(fxmatpar(iseg2).GT.0.d0) THEN
              fp23=0.d0
              IF(fp32.GT. fxmatpar(iseg2)) fp32=fxmatpar(iseg2)
            ELSE
              fp32=0.d0
              IF(fp23.GT.-fxmatpar(iseg2)) fp23=-fxmatpar(iseg2)
            ENDIF
          ELSE
            IF(fxmatpar(iseg2).GT.0.d0) THEN
              fp32=0.d0
              IF(fp23.GT. fxmatpar(iseg2)) fp23=fxmatpar(iseg2)
            ELSE
              fp23=0.d0
              IF(fp32.GT.-fxmatpar(iseg2)) fp32=-fxmatpar(iseg2)
            ENDIF
          ENDIF
!
!         SEGMENT 3
!
          IF(oriseg(ielem,3).EQ.1) THEN
            IF(fxmatpar(iseg3).GT.0.d0) THEN
              fp31=0.d0
              IF(fp13.GT. fxmatpar(iseg3)) fp13=fxmatpar(iseg3)
            ELSE
              fp13=0.d0
              IF(fp31.GT.-fxmatpar(iseg3)) fp31=-fxmatpar(iseg3)
            ENDIF
          ELSE
            IF(fxmatpar(iseg3).GT.0.d0) THEN
              fp13=0.d0
              IF(fp31.GT.fxmatpar(iseg3)) fp31=fxmatpar(iseg3)
            ELSE
              fp31=0.d0
              IF(fp13.GT.-fxmatpar(iseg3)) fp13=-fxmatpar(iseg3)
            ENDIF
          ENDIF
!
!         END OF CORRECTION OF FLUXES
!
          fn1=fn%R(i1)
          fn2=fn%R(i2)
          fn3=fn%R(i3)
!
!         MINIMUM OF THE 1-TETA
!
          mt1=1.d0-teta(i1)
          mt2=1.d0-teta(i2)
          mt3=1.d0-teta(i3)
          min12=min(mt1,mt2)
          min13=min(mt1,mt3)
          min23=min(mt2,mt3)
!
!         PART OF CONTRIBUTIONS THAT WILL NOT BE LIMITED, IMMEDIATELY ASSEMBLED
!
          fi_i(i1)=fi_i(i1)+fp12*(fn1*(mt1-min12)-fn2*(mt2-min12))
     &                     +fp13*(fn1*(mt1-min13)-fn3*(mt3-min13))
          fi_i(i2)=fi_i(i2)+fp21*(fn2*(mt2-min12)-fn1*(mt1-min12))
     &                     +fp23*(fn2*(mt2-min23)-fn3*(mt3-min23))
          fi_i(i3)=fi_i(i3)+fp31*(fn3*(mt3-min13)-fn1*(mt1-min13))
     &                     +fp32*(fn3*(mt3-min23)-fn2*(mt2-min23))
!
!         NOW PART OF CONTRIBUTIONS THAT WILL BE LIMITED
!
          coef=su(ielem)*tiers
!
!         AS CLASSICAL N SCHEME, BUT WITH DERIVATIVE IN TIME ADDED
!         HDFDT MUST BE ((1.D0-TETA)*H+TETA*HN)*(FSTAR-F)/DDT
!
          fincorr1=
     &    hdfdt(i1)*coef+fp12*(fn1-fn2)*min12+fp13*(fn1-fn3)*min13
          fincorr2=
     &    hdfdt(i2)*coef+fp21*(fn2-fn1)*min12+fp23*(fn2-fn3)*min23
          fincorr3=
     &    hdfdt(i3)*coef+fp31*(fn3-fn1)*min13+fp32*(fn3-fn2)*min23
!
!         PHITCOR IS THE NEW TOTAL RESIDUAL,
!
          phitcor=fincorr1+fincorr2+fincorr3
!
          IF(phitcor.GT.epsphi) THEN
!           PSI REDUCTION
            beta1=max(fincorr1,0.d0)
            beta2=max(fincorr2,0.d0)
            beta3=max(fincorr3,0.d0)
            coef=phitcor/(beta1+beta2+beta3)
            fi_i(i1)=fi_i(i1)+beta1*coef
            fi_i(i2)=fi_i(i2)+beta2*coef
            fi_i(i3)=fi_i(i3)+beta3*coef
          ELSEIF(phitcor.LT.-epsphi) THEN
!           PSI REDUCTION
            beta1=min(fincorr1,0.d0)
            beta2=min(fincorr2,0.d0)
            beta3=min(fincorr3,0.d0)
            coef=phitcor/(beta1+beta2+beta3)
            fi_i(i1)=fi_i(i1)+beta1*coef
            fi_i(i2)=fi_i(i2)+beta2*coef
            fi_i(i3)=fi_i(i3)+beta3*coef
          ELSE
!           NO REDUCTION
            fi_i(i1)=fi_i(i1)+fincorr1
            fi_i(i2)=fi_i(i2)+fincorr2
            fi_i(i3)=fi_i(i3)+fincorr3
          ENDIF
!
        ENDDO
!
      ELSEIF(yaflulim) THEN
!
!       THE SAME BUT WITH FLUX LIMITATION GIVEN
!
        DO ielem = 1,nelem
!
          k1 = -phiel(ielem,1)
          k2 = -phiel(ielem,2)
          k3 = -phiel(ielem,3)
!
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
!
          iseg1=eltseg(ielem,1)
          iseg2=eltseg(ielem,2)
          iseg3=eltseg(ielem,3)
!
!         STARTS WITH N-SCHEME (EQUIVALENT TO LEO POSTMA'S IMPLEMENTATION)
!         FIJ HERE LIKE LAMBDA(I,J) IN BOOK
!         I.E. FIJ>0 MEANS FLUX FROM J TO I !!!!!
!
          fp12=max(min(k1,-k2),0.d0)*flulim(iseg1)
          fp23=max(min(k2,-k3),0.d0)*flulim(iseg2)
          fp31=max(min(k3,-k1),0.d0)*flulim(iseg3)
          fp21=max(min(k2,-k1),0.d0)*flulim(iseg1)
          fp32=max(min(k3,-k2),0.d0)*flulim(iseg2)
          fp13=max(min(k1,-k3),0.d0)*flulim(iseg3)
!
!         CORRECTING THE FLUXES WHEN THEIR SIGN IS NOT THE SAME
!         AS THE ASSEMBLED VALUE
!
!         SEGMENT 1
!
          IF(oriseg(ielem,1).EQ.1) THEN
            IF(fxmatpar(iseg1).GT.0.d0) THEN
              fp12=0.d0
              IF(fp21.GT. fxmatpar(iseg1)) fp21=fxmatpar(iseg1)
            ELSE
              fp21=0.d0
              IF(fp12.GT.-fxmatpar(iseg1)) fp12=-fxmatpar(iseg1)
            ENDIF
          ELSE
            IF(fxmatpar(iseg1).GT.0.d0) THEN
              fp21=0.d0
              IF(fp12.GT. fxmatpar(iseg1)) fp12=fxmatpar(iseg1)
            ELSE
              fp12=0.d0
              IF(fp21.GT.-fxmatpar(iseg1)) fp21=-fxmatpar(iseg1)
            ENDIF
          ENDIF
!
!         SEGMENT 2
!
          IF(oriseg(ielem,2).EQ.1) THEN
            IF(fxmatpar(iseg2).GT.0.d0) THEN
              fp23=0.d0
              IF(fp32.GT. fxmatpar(iseg2)) fp32=fxmatpar(iseg2)
            ELSE
              fp32=0.d0
              IF(fp23.GT.-fxmatpar(iseg2)) fp23=-fxmatpar(iseg2)
            ENDIF
          ELSE
            IF(fxmatpar(iseg2).GT.0.d0) THEN
              fp32=0.d0
              IF(fp23.GT. fxmatpar(iseg2)) fp23=fxmatpar(iseg2)
            ELSE
              fp23=0.d0
              IF(fp32.GT.-fxmatpar(iseg2)) fp32=-fxmatpar(iseg2)
            ENDIF
          ENDIF
!
!         SEGMENT 3
!
          IF(oriseg(ielem,3).EQ.1) THEN
            IF(fxmatpar(iseg3).GT.0.d0) THEN
              fp31=0.d0
              IF(fp13.GT. fxmatpar(iseg3)) fp13=fxmatpar(iseg3)
            ELSE
              fp13=0.d0
              IF(fp31.GT.-fxmatpar(iseg3)) fp31=-fxmatpar(iseg3)
            ENDIF
          ELSE
            IF(fxmatpar(iseg3).GT.0.d0) THEN
              fp13=0.d0
              IF(fp31.GT.fxmatpar(iseg3)) fp31=fxmatpar(iseg3)
            ELSE
              fp31=0.d0
              IF(fp13.GT.-fxmatpar(iseg3)) fp13=-fxmatpar(iseg3)
            ENDIF
          ENDIF
!
!         END OF CORRECTION OF FLUXES
!
          fn1=fn%R(i1)
          fn2=fn%R(i2)
          fn3=fn%R(i3)
!
!         MINIMUM OF THE 1-TETA
!
          mt1=1.d0-teta(i1)
          mt2=1.d0-teta(i2)
          mt3=1.d0-teta(i3)
          min12=min(mt1,mt2)
          min13=min(mt1,mt3)
          min23=min(mt2,mt3)
!
!         PART OF CONTRIBUTIONS THAT WILL NOT BE LIMITED, IMMEDIATELY ASSEMBLED
!
          fi_i(i1)=fi_i(i1)+fp12*(fn1*(mt1-min12)-fn2*(mt2-min12))
     &                     +fp13*(fn1*(mt1-min13)-fn3*(mt3-min13))
          fi_i(i2)=fi_i(i2)+fp21*(fn2*(mt2-min12)-fn1*(mt1-min12))
     &                     +fp23*(fn2*(mt2-min23)-fn3*(mt3-min23))
          fi_i(i3)=fi_i(i3)+fp31*(fn3*(mt3-min13)-fn1*(mt1-min13))
     &                     +fp32*(fn3*(mt3-min23)-fn2*(mt2-min23))
!
!         NOW PART OF CONTRIBUTIONS THAT WILL BE LIMITED
!
          coef=su(ielem)*tiers
!
!         AS CLASSICAL N SCHEME, BUT WITH DERIVATIVE IN TIME ADDED
!         HDFDT MUST BE ((1.D0-TETA)*H+TETA*HN)*(FSTAR-F)/DDT
!
          fincorr1=hdfdt(i1)*coef
     &            +fp12*(fn1-fn2)*min12+fp13*(fn1-fn3)*min13
          fincorr2=hdfdt(i2)*coef
     &            +fp21*(fn2-fn1)*min12+fp23*(fn2-fn3)*min23
          fincorr3=hdfdt(i3)*coef
     &            +fp31*(fn3-fn1)*min13+fp32*(fn3-fn2)*min23
!
!         PHITCOR IS THE NEW TOTAL RESIDUAL,
!
          phitcor=fincorr1+fincorr2+fincorr3
!
          IF(phitcor.GT.epsphi) THEN
!           PSI REDUCTION
            beta1=max(fincorr1,0.d0)
            beta2=max(fincorr2,0.d0)
            beta3=max(fincorr3,0.d0)
            coef=phitcor/(beta1+beta2+beta3)
            fi_i(i1)=fi_i(i1)+beta1*coef
            fi_i(i2)=fi_i(i2)+beta2*coef
            fi_i(i3)=fi_i(i3)+beta3*coef
          ELSEIF(phitcor.LT.-epsphi) THEN
!           PSI REDUCTION
            beta1=min(fincorr1,0.d0)
            beta2=min(fincorr2,0.d0)
            beta3=min(fincorr3,0.d0)
            coef=phitcor/(beta1+beta2+beta3)
            fi_i(i1)=fi_i(i1)+beta1*coef
            fi_i(i2)=fi_i(i2)+beta2*coef
            fi_i(i3)=fi_i(i3)+beta3*coef
          ELSE
!           NO REDUCTION
            fi_i(i1)=fi_i(i1)+fincorr1
            fi_i(i2)=fi_i(i2)+fincorr2
            fi_i(i3)=fi_i(i3)+fincorr3
          ENDIF
!
        ENDDO
!
      ELSE
!
!       HERE NO FLUX LIMITATION
!
        DO ielem = 1,nelem
!
          k1 = -phiel(ielem,1)
          k2 = -phiel(ielem,2)
          k3 = -phiel(ielem,3)
!
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
!
          iseg1=eltseg(ielem,1)
          iseg2=eltseg(ielem,2)
          iseg3=eltseg(ielem,3)
!
!         STARTS WITH N-SCHEME (EQUIVALENT TO LEO POSTMA'S IMPLEMENTATION)
!         FIJ HERE LIKE LAMBDA(I,J) IN BOOK
!         I.E. FIJ>0 MEANS FLUX FROM J TO I !!!!!
!
          fp12=max(min(k1,-k2),0.d0)
          fp23=max(min(k2,-k3),0.d0)
          fp31=max(min(k3,-k1),0.d0)
          fp21=max(min(k2,-k1),0.d0)
          fp32=max(min(k3,-k2),0.d0)
          fp13=max(min(k1,-k3),0.d0)
!
!         CORRECTING THE FLUXES WHEN THEIR SIGN IS NOT THE SAME
!         AS THE ASSEMBLED VALUE, KNOWING THAT ALL THE FPIJ ARE
!         POSITIVE BY CONSTRUCTION
!
!         SEGMENT 1
!
          IF(oriseg(ielem,1).EQ.1) THEN
            IF(fxmatpar(iseg1).GT.0.d0) THEN
              fp12=0.d0
              IF(fp21.GT. fxmatpar(iseg1)) fp21=fxmatpar(iseg1)
            ELSE
              fp21=0.d0
              IF(fp12.GT.-fxmatpar(iseg1)) fp12=-fxmatpar(iseg1)
            ENDIF
          ELSE
            IF(fxmatpar(iseg1).GT.0.d0) THEN
              fp21=0.d0
              IF(fp12.GT. fxmatpar(iseg1)) fp12=fxmatpar(iseg1)
            ELSE
              fp12=0.d0
              IF(fp21.GT.-fxmatpar(iseg1)) fp21=-fxmatpar(iseg1)
            ENDIF
          ENDIF
!
!         SEGMENT 2
!
          IF(oriseg(ielem,2).EQ.1) THEN
            IF(fxmatpar(iseg2).GT.0.d0) THEN
              fp23=0.d0
              IF(fp32.GT. fxmatpar(iseg2)) fp32=fxmatpar(iseg2)
            ELSE
              fp32=0.d0
              IF(fp23.GT.-fxmatpar(iseg2)) fp23=-fxmatpar(iseg2)
            ENDIF
          ELSE
            IF(fxmatpar(iseg2).GT.0.d0) THEN
              fp32=0.d0
              IF(fp23.GT. fxmatpar(iseg2)) fp23=fxmatpar(iseg2)
            ELSE
              fp23=0.d0
              IF(fp32.GT.-fxmatpar(iseg2)) fp32=-fxmatpar(iseg2)
            ENDIF
          ENDIF
!
!         SEGMENT 3
!
          IF(oriseg(ielem,3).EQ.1) THEN
            IF(fxmatpar(iseg3).GT.0.d0) THEN
              fp31=0.d0
              IF(fp13.GT. fxmatpar(iseg3)) fp13=fxmatpar(iseg3)
            ELSE
              fp13=0.d0
              IF(fp31.GT.-fxmatpar(iseg3)) fp31=-fxmatpar(iseg3)
            ENDIF
          ELSE
            IF(fxmatpar(iseg3).GT.0.d0) THEN
              fp13=0.d0
              IF(fp31.GT.fxmatpar(iseg3)) fp31=fxmatpar(iseg3)
            ELSE
              fp31=0.d0
              IF(fp13.GT.-fxmatpar(iseg3)) fp13=-fxmatpar(iseg3)
            ENDIF
          ENDIF
!
!         END OF CORRECTION OF FLUXES
!
          fn1=fn%R(i1)
          fn2=fn%R(i2)
          fn3=fn%R(i3)
!
!         MINIMUM OF THE 1-TETA
!
          mt1=1.d0-teta(i1)
          mt2=1.d0-teta(i2)
          mt3=1.d0-teta(i3)
          min12=min(mt1,mt2)
          min13=min(mt1,mt3)
          min23=min(mt2,mt3)
!
!         PART OF CONTRIBUTIONS THAT WILL NOT BE LIMITED, IMMEDIATELY ASSEMBLED
!
          fi_i(i1)=fi_i(i1)+fp12*(fn1*(mt1-min12)-fn2*(mt2-min12))
     &                     +fp13*(fn1*(mt1-min13)-fn3*(mt3-min13))
          fi_i(i2)=fi_i(i2)+fp21*(fn2*(mt2-min12)-fn1*(mt1-min12))
     &                     +fp23*(fn2*(mt2-min23)-fn3*(mt3-min23))
          fi_i(i3)=fi_i(i3)+fp31*(fn3*(mt3-min13)-fn1*(mt1-min13))
     &                     +fp32*(fn3*(mt3-min23)-fn2*(mt2-min23))
!
!         NOW PART OF CONTRIBUTIONS THAT WILL BE LIMITED
!
          coef=su(ielem)*tiers
!
!         AS CLASSICAL N SCHEME, BUT WITH DERIVATIVE IN TIME ADDED
!         HDFDT MUST BE ((1.D0-TETA)*H+TETA*HN)*(FSTAR-F)/DDT
!
          fincorr1=
     &    hdfdt(i1)*coef+fp12*(fn1-fn2)*min12+fp13*(fn1-fn3)*min13
          fincorr2=
     &    hdfdt(i2)*coef+fp21*(fn2-fn1)*min12+fp23*(fn2-fn3)*min23
          fincorr3=
     &    hdfdt(i3)*coef+fp31*(fn3-fn1)*min13+fp32*(fn3-fn2)*min23
!
!         PHITCOR IS THE NEW TOTAL RESIDUAL,
!
          phitcor=fincorr1+fincorr2+fincorr3
!
          IF(phitcor.GT.epsphi) THEN
!           PSI REDUCTION
            beta1=max(fincorr1,0.d0)
            beta2=max(fincorr2,0.d0)
            beta3=max(fincorr3,0.d0)
            coef=phitcor/(beta1+beta2+beta3)
            fi_i(i1)=fi_i(i1)+beta1*coef
            fi_i(i2)=fi_i(i2)+beta2*coef
            fi_i(i3)=fi_i(i3)+beta3*coef
          ELSEIF(phitcor.LT.-epsphi) THEN
!           PSI REDUCTION
            beta1=min(fincorr1,0.d0)
            beta2=min(fincorr2,0.d0)
            beta3=min(fincorr3,0.d0)
            coef=phitcor/(beta1+beta2+beta3)
            fi_i(i1)=fi_i(i1)+beta1*coef
            fi_i(i2)=fi_i(i2)+beta2*coef
            fi_i(i3)=fi_i(i3)+beta3*coef
          ELSE
!           NO REDUCTION
            fi_i(i1)=fi_i(i1)+fincorr1
            fi_i(i2)=fi_i(i2)+fincorr2
            fi_i(i3)=fi_i(i3)+fincorr3
          ENDIF
!
        ENDDO
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END