tvf_imp.F

Go to the documentation of this file.

!                   ******************
                    SUBROUTINE tvf_imp
!                   ******************
!
     &(f,fc,fxmat,fxmatpar,
     & dt,fxbor,hnp1mt,fbor,smh,yasmh,fscexp,
     & nseg,npoin,nptfr,gloseg,sizglo,nbor,limtra,kdir,kddl,optsou,
     & iopt2,flbortra,surnit,mesh,sf,rain,pluie,train,tetaf,infogt,
     & volu2d,sm,psiexp,am2,bb,slvpsi,
     & predictor,corrector,icor,ncor,massou)
!
!***********************************************************************
! BIEF   V7P2
!***********************************************************************
!
!brief    Semi-implicit distributive scheme.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        26/06/2015
!+        V7P1
!+   First version.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        01/09/2015
!+        V7P1
!+   Optimisation. Matrix simplified into a diagonal except for the last
!+   correction. Does not spoil mass conservation and results virtually
!+   unchanged.
!+   Jacobi solver programmed here, when GMRES or direct solver
!+   is not asked.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        30/06/2016
!+        V7P2
!+   Computing the global extrema for clipping the results that may be
!+   out of range due to truncation errors. True errors may be hidden
!+   but will be seen with the mass balance.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        30/06/2016
!+        V7P2
!+   Relaxing the clipping for monotony in case of rain (evaporation
!+   may break the monotony).
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DT             |-->| TIME-STEP
!| F              |<--| VALUES OF F AT TIME N+1 OF SUB-ITERATION
!| FBOR           |-->| VALUES OF F AT THE PRESCRIBED BOUNDARIES
!| FC             |-->| VALUES OF F AT TIME N OF SUB-ITERATION
!| FLBORTRA       |<->| FLUX OF TRACER AT THE BOUNDARIES
!| FSCEXP         |-->| EXPLICIT SOURCE TERM FOR F
!| FXBOR          |-->| FLUXES ON BOUNDARIES
!| FXMAT          |-->| FLUXES (NON ASSEMBLED IN PARALLEL)
!| FXMATPAR       |-->| FLUXES (ASSEMBLED IN PARALLEL)
!| GLOSEG         |-->| GLOBAL NUMBER OF THE 2 POINTS OF A SEGMENT.
!| ICOR           |-->| CURRENT CORRECTION NUMBER
!| INFOGT         |-->| IF YES, PRINT INFORMATION ON SOLVER
!| IOPT2          |-->| 0: CONSERVATIVE ADVECTION FIELD
!|                |   | 1: NON CONSERVATIVE ADVECTION FIELD
!| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
!| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
!| LIMTRA         |-->| TECHNICAL BOUNDARY CONDITIONS FOR TRACERS
!| MASSOU         |-->| MASS ADDED BY SOURCE TERM
!| MESH           |-->| MESH STRUCTURE
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| NCOR           |-->| TOTAL NUMBER OF CORRECTIONS ASKED
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NSEG           |-->| NUMBER OF SEGMENTS
!| OPTSOU         |-->| TYPE OF SOURCES
!|                |   | 1: NORMAL
!|                |   | 2: DIRAC
!| PLUIE          |-->| RAIN OR EVAPORATION, IN M/S
!| RAIN           |-->| IF YES: RAIN OR EVAPORATION
!| SF             |<->| BIEF_OBJ STRUCTURE OF F
!| SIZGLO         |-->| FIRST DIMENSION OF GLOSEG
!| SMH            |-->| SOURCE TERM IN CONTINUITY EQUATION
!| SURNIT         |-->| SURNIT=1/NIT
!| TETAF          |-->| IMPLICITATION COEFFICIENT ON F
!| HNP1MT         |-->| INTERMEDIATE DEPTH H(N+1-TETA)
!| TRAIN          |-->| VALUE OF TRACER IN RAIN
!| VOLU2D         |-->| INTEGRALS OF TEST FUNCTIONS, NOT ASSEMBLED IN //
!| YASMH          |-->| IF YES, SMH MUST BE TAKEN INTO ACCOUNT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_tvf_imp => tvf_imp
      USE declarations_special
      USE interface_parallel
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NSEG,NPOIN,NPTFR,KDIR,KDDL
      INTEGER, INTENT(IN)             :: SIZGLO,OPTSOU,IOPT2
      INTEGER, INTENT(IN)             :: ICOR,NCOR
      INTEGER, INTENT(IN)             :: GLOSEG(sizglo,2)
      INTEGER, INTENT(IN)             :: NBOR(nptfr),LIMTRA(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DT,SURNIT,TRAIN,TETAF(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: FLBORTRA(nptfr),MASSOU
      DOUBLE PRECISION, INTENT(INOUT) :: F(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: FXBOR(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: FC(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: HNP1MT(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: SMH(*)
      DOUBLE PRECISION, INTENT(IN)    :: VOLU2D(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: PLUIE(*),PSIEXP(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: FSCEXP(*)
      DOUBLE PRECISION, INTENT(IN)    :: FBOR(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: FXMAT(nseg),FXMATPAR(nseg)
      LOGICAL, INTENT(IN)             :: YASMH,RAIN,INFOGT
      LOGICAL, INTENT(IN)             :: PREDICTOR,CORRECTOR
      TYPE(bief_obj),  INTENT(INOUT)  :: SF,SM,AM2,BB
      TYPE(bief_mesh), INTENT(INOUT)  :: MESH
      TYPE(slvcfg),    INTENT(INOUT)  :: SLVPSI
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,N,I1,I2
      DOUBLE PRECISION NORMR,NORMB,FMIN,FMAX
      TYPE(bief_obj), POINTER :: BB1,SURDIAG,R
!
      INTRINSIC sqrt
!
!-----------------------------------------------------------------------
!
      bb1    =>bb%ADR(1)%P
      r      =>bb%ADR(2)%P
      surdiag=>bb%ADR(3)%P
      CALL cpstvc(sf,bb1)
      CALL cpstvc(sf,r)
      CALL cpstvc(sf,surdiag)
!
!-----------------------------------------------------------------------
!
!     COMPUTING GLOBAL EXTREMA, FOR CLIPPING TRUNCATION ERRORS
!
      fmin=fc(1)
      fmax=fc(1)
      DO i=1,npoin
        fmin=min(fmin,fc(i))
        fmax=max(fmax,fc(i))
      ENDDO
      DO i=1,nptfr
        IF(limtra(i).EQ.kdir) THEN
          fmin=min(fmin,fbor(i))
          fmax=max(fmax,fbor(i))
        ENDIF
      ENDDO
      IF(yasmh) THEN
        DO i=1,npoin
          fmin=min(fmin,fscexp(i))
          fmax=max(fmax,fscexp(i))
        ENDDO
      ENDIF
      IF(rain) THEN
        fmin=min(fmin,train)
        fmax=max(fmax,train)
      ENDIF
      IF(ncsize.GT.1) THEN
        fmin=p_min(fmin)
        fmax=p_max(fmax)
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(predictor) THEN
!
        IF(iopt2.EQ.0) THEN
!         CONSERVATIVE ADVECTION FIELD
          DO i = 1,npoin
            f(i) = fc(i)
          ENDDO
        ELSE
          WRITE(lu,*) 'TVF_IMP: UNKNOWN OPTION'
          CALL plante(1)
          stop
        ENDIF
!
      ENDIF
!
!     BUILDING THE MATRIX AND THE RIGHT-HAND SIDE
!     THE MATRIX IS NOT DONE WITH A CALL MATRIX
!     SO ITS FEATURES HAVE TO BE HARDCODED HERE
!
      am2%X%DIM1=2*nseg
      am2%X%DIM2=1
      am2%TYPEXT='Q'
      am2%TYPDIA='Q'
      am2%ELMLIN=11
      am2%ELMCOL=11
      CALL cpstvc(sf,am2%D)
!
!     DIAGONAL OF MATRIX
!
      DO i=1,npoin
        am2%D%R(i)=hnp1mt(i)*volu2d(i)
      ENDDO
!
!     RIGHT HAND SIDE
!
!     TERM FROM THE DERIVATIVE IN TIME
      IF(predictor) THEN
        DO i=1,npoin
          sm%R(i)=am2%D%R(i)*fc(i)
        ENDDO
      ELSEIF(corrector) THEN
!       THE PREDICTOR IS TAKEN, AT THIS LEVEL IT IS STILL F
        DO i=1,npoin
          sm%R(i)=am2%D%R(i)*f(i)
        ENDDO
      ELSE
        WRITE(lu,*) 'TVF_IMP, CHECK ARGUMENTS PREDICTOR, CORRECTOR'
        CALL plante(1)
        stop
      ENDIF
!
!     IMPLICIT BOUNDARY TERM
!
      DO i=1,nptfr
        n=nbor(i)
        IF(limtra(i).EQ.kdir) THEN
          am2%D%R(n)=am2%D%R(n)-dt*tetaf(n)*fxbor(i)
        ENDIF
      ENDDO
!
!     DIAGONAL AND OFF-DIAGONAL TERMS
!
      IF(icor.LT.ncor) THEN
!       SYSTEM SIMPLIFIED, MASS SPOILED, BUT IT WILL BE CORRECTED AFTER
!       BY THE LAST CORRECTION
!       HERE EXTRA-DIAGONAL TERMS NOT BUILT
        DO i=1,nseg
          i1=gloseg(i,1)
          i2=gloseg(i,2)
          IF(fxmatpar(i).LT.0.d0) THEN
            am2%D%R(i1) = am2%D%R(i1) - dt*tetaf(i1)*fxmat(i)
          ELSE
            am2%D%R(i2) = am2%D%R(i2) + dt*tetaf(i2)*fxmat(i)
          ENDIF
        ENDDO
      ELSE
!       NO SIMPLIFICATION, REAL MATRIX
        DO i=1,nseg
          i1=gloseg(i,1)
          i2=gloseg(i,2)
          IF(fxmatpar(i).LT.0.d0) THEN
            am2%D%R(i1) = am2%D%R(i1) - dt*tetaf(i1)*fxmat(i)
            am2%X%R(i)=dt*tetaf(i2)*fxmat(i)
            am2%X%R(i+nseg)=0.d0
          ELSE
            am2%D%R(i2) = am2%D%R(i2) + dt*tetaf(i2)*fxmat(i)
            am2%X%R(i)=0.d0
            am2%X%R(i+nseg)=-dt*tetaf(i1)*fxmat(i)
          ENDIF
        ENDDO
!
      ENDIF
!
!     SOURCES IN CONTINUITY EQUATION (SMH)
!
      IF(yasmh) THEN
        IF(optsou.EQ.1) THEN
!         SEE DIFSOU, FSCEXP CONTAINS THE VALUE OF THE TRACER AT THE SOURCE
          DO i=1,npoin
            am2%D%R(i)=am2%D%R(i)
     &                +dt*tetaf(i)*volu2d(i)*max(smh(i),0.d0)
            sm%R(i)=sm%R(i)+(fscexp(i)-(1.d0-tetaf(i))*fc(i))
     &                     *volu2d(i)*dt*max(smh(i),0.d0)
          ENDDO
        ELSEIF(optsou.EQ.2) THEN
          DO i=1,npoin
            am2%D%R(i)=am2%D%R(i)+dt*tetaf(i)*max(smh(i),0.d0)
            sm%R(i)=sm%R(i)+(fscexp(i)-(1.d0-tetaf(i))*fc(i))
     &                     *dt*max(smh(i),0.d0)
          ENDDO
        ENDIF
      ENDIF
!
!     TERMES BII * CIN ET BIJ * CJN
!
!     ASSEMBLED CONTRIBUTION, DONE BY FLUX_EF_VF_3, POSSIBLY
!     WITH DERIVATIVE IN TIME.
      DO i=1,npoin
        sm%R(i)=sm%R(i)-dt*psiexp(i)
      ENDDO
!
!     ADD FLUX ON BOUNDARY AND OTHER TERMS
!
      DO i=1,nptfr
        IF(limtra(i).EQ.kdir) THEN
          n=nbor(i)
          sm%R(n)=sm%R(n)+dt*fxbor(i)*((1.d0-tetaf(n))*fc(n)-fbor(i))
        ENDIF
      ENDDO
!
!     RAIN-EVAPORATION
!
      IF(rain) THEN
        DO i=1,npoin
          IF(pluie(i).GT.0.d0) THEN
!           REAL RAIN, VALUE IN RAIN CONSIDERED...
            sm%R(i)=sm%R(i)+dt*volu2d(i)*pluie(i)
     &          *(train-(1.d0-tetaf(i))*fc(i))
          ELSE
!           EVAPORATION, VALUE IN RAIN NOT CONSIDERED...
            sm%R(i)=sm%R(i)+dt*volu2d(i)*pluie(i)
     &          *(     -(1.d0-tetaf(i))*fc(i))
          ENDIF
          am2%D%R(i)=am2%D%R(i)+dt*tetaf(i)*volu2d(i)*pluie(i)
        ENDDO
      ENDIF
!
!     TIDAL FLATS
!
      CALL os('X=Y     ',x=bb1,y=am2%D)
      IF(ncsize.GT.1) CALL parcom(bb1,2,mesh)
      DO i=1,npoin
!       SEE PRECD1 EPSILON HERE MUST BE GREATER THAN 1.D-20, OR PRECD1 WILL
!       DO THE CLIPPING ITSELF, IN A LESS CONSISTANT WAY
!       THE TEST MUST BE DONE ON THE ASSEMBLED DIAGONAL
        IF(bb1%R(i).LT.1.d-15) THEN
!         VOLU2D IS A NON ASSEMBLED COEFFICIENT, ANY SUCH COEFFICIENT WOULD WORK...
!         DRY POINT THAT RECEIVES NO WATER, F=FC IS GIVEN AS EQUATION
          am2%D%R(i)=volu2d(i)
          sm%R(i)   =volu2d(i)*fc(i)
        ENDIF
      ENDDO
!
!     SOLVING THE FINAL LINEAR SYSTEM
!
      IF(icor.LT.ncor) THEN
!       HERE ONE ITERATION OF JACOBI, OFF-DIAGONAL TERMS BUILT ON THE SPOT
        DO i=1,nseg
          i1=gloseg(i,1)
          i2=gloseg(i,2)
          IF(fxmatpar(i).LT.0.d0) THEN
            sm%R(i1)=sm%R(i1)-dt*tetaf(i2)*fxmat(i)*f(i2)
          ELSE
            sm%R(i2)=sm%R(i2)+dt*tetaf(i1)*fxmat(i)*f(i1)
          ENDIF
        ENDDO
        IF(ncsize.GT.1) THEN
          CALL parcom(am2%D,2,mesh)
          CALL parcom(sm,2,mesh)
        ENDIF
        DO i=1,npoin
          f(i)=sm%R(i)/am2%D%R(i)
        ENDDO
      ELSE
        IF(slvpsi%SLV.EQ.7.OR.slvpsi%SLV.EQ.8) THEN
#if defined COMPAD
          CALL ad_solve(sf,am2,sm,bb,slvpsi,infogt,mesh,am2)
#else
          CALL solve(sf,am2,sm,bb,slvpsi,infogt,mesh,am2)
#endif
        ELSE
!         IF NOT GMRES OR DIRECT, JACOBI SOLUTION
          n=0
          CALL os('X=Y     ',x=bb1,y=sm)
          CALL os('X=Y     ',x=surdiag,y=am2%D)
          IF(ncsize.GT.1) THEN
            CALL parcom(surdiag,2,mesh)
            CALL parcom(bb1,2,mesh)
          ENDIF
          CALL os('X=1/Y   ',x=surdiag,y=surdiag)
          normb=sqrt(p_dots(bb1,bb1,mesh))
100       CONTINUE
          n=n+1
!         ONE ITERATION OF JACOBI
          CALL os('X=Y     ',x=bb1,y=sm)
          DO i=1,nseg
            i1=gloseg(i,1)
            i2=gloseg(i,2)
            IF(fxmatpar(i).LT.0.d0) THEN
              bb1%R(i1)=bb1%R(i1)-am2%X%R(i     )*f(i2)
            ELSE
              bb1%R(i2)=bb1%R(i2)-am2%X%R(i+nseg)*f(i1)
            ENDIF
          ENDDO
          IF(ncsize.GT.1) CALL parcom(bb1,2,mesh)
!         NEW SOLUTION IN BB1
          DO i=1,npoin
            bb1%R(i)=bb1%R(i)*surdiag%R(i)
          ENDDO
!         END OF ONE ITERATION, COMPUTING THE RESIDUAL
!         AX-B = (A%D+A%X)X-B = A%X*(NEW F - OLD F)
          CALL os('X=0     ',x=r)
          DO i=1,nseg
            i1=gloseg(i,1)
            i2=gloseg(i,2)
            IF(fxmatpar(i).LT.0.d0) THEN
              r%R(i1)=r%R(i1)-am2%X%R(i     )*(f(i2)-bb1%R(i2))
            ELSE
              r%R(i2)=r%R(i2)-am2%X%R(i+nseg)*(f(i1)-bb1%R(i1))
            ENDIF
          ENDDO
          IF(ncsize.GT.1) CALL parcom(r,2,mesh)
          normr=sqrt(p_dots(r,r,mesh))
!         COPY OF NEW SOLUTION ON F, WITH CLIPPING WITH GLOBAL EXTREMA
!         TO COPE WITH TRUNCATION ERRORS. IF CLIPPING TRUE ERRORS IT
!         WILL DO MASS ERRORS
          IF(rain) THEN
!           NO CLIPPING HERE SINCE EVAPORATION MAY BREAK MONOTONY
            DO i=1,npoin
              f(i)=bb1%R(i)
            ENDDO
          ELSE
            DO i=1,npoin
              f(i)=max(min(bb1%R(i),fmax),fmin)
            ENDDO
          ENDIF
          IF(n.LT.slvpsi%NITMAX.AND.normr.GT.slvpsi%EPS*normb) THEN
            GO TO 100
          ENDIF
!         DONE !
          IF(infogt) THEN
            IF(normr.GT.slvpsi%EPS*normb) THEN
              IF(normb.GT.1.d0) THEN
                WRITE(lu,*)
     &            'JACOBI: MAXIMUM ITERATIONS REACHED ',n,
     &            ' ITERATIONS, RELATIVE PRECISION = ',normr/normb
              ELSE
                WRITE(lu,*)
     &            'JACOBI: MAXIMUM ITERATIONS REACHED ',n,
     &            ' ITERATIONS, ABSOLUTE PRECISION = ',normr
              ENDIF
            ELSE
              IF(normb.GT.1.d0) THEN
                WRITE(lu,*) 'JACOBI:',n,
     &            ' ITERATIONS, RELATIVE PRECISION = ',normr/normb
              ELSE
                WRITE(lu,*) 'JACOBI:',n,
     &            ' ITERATIONS, ABSOLUTE PRECISION = ',normr
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
!
!     ON EXITS, EXITING FLUX DEPENDS ON F AND FN AT EVERY SUB-TIME
!     STEP, SO IT MUST BE COMPUTED AT THIS LEVEL.
!     THE CASE KDIR IS TREATED IN CVTRVF
!
!     FLUX AND ADDED MASS FOR MASS BALANCE
!
      IF(icor.EQ.ncor) THEN
!       BOUNDARIES
        DO i=1,nptfr
          IF(limtra(i).EQ.kdir) THEN
            flbortra(i)=flbortra(i)+fxbor(i)*fbor(i)*surnit
          ELSEIF(limtra(i).EQ.kddl) THEN
            n=nbor(i)
            flbortra(i)=flbortra(i)
     &            +fxbor(i)*(tetaf(n)*f(n)+(1.d0-tetaf(n))*fc(n))*surnit
          ENDIF
        ENDDO
!       SOURCES
        IF(yasmh) THEN
!         FOR MASS BALANCE
          IF(optsou.EQ.1) THEN
            DO i=1,mesh%NPOIN
              IF(smh(i).GE.0.d0) THEN
                massou=massou+volu2d(i)*fscexp(i)*dt*smh(i)
              ELSE
                massou=massou+volu2d(i)*
     &                (tetaf(i)*f(i)+(1.d0-tetaf(i))*fc(i))*dt*smh(i)
              ENDIF
            ENDDO
          ELSEIF(optsou.EQ.2) THEN
            DO i=1,mesh%NPOIN
              IF(smh(i).GE.0.d0) THEN
                massou=massou+fscexp(i)*dt*smh(i)
              ELSE
                massou=massou+
     &                  (tetaf(i)*f(i)+(1.d0-tetaf(i))*fc(i))*dt*smh(i)
              ENDIF
            ENDDO
          ENDIF
        ENDIF
!       RAIN
!       DONE IN BILANT FOR ALL SCHEMES
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END