tvf.f

Go to the documentation of this file.

!                   **************
                    SUBROUTINE tvf
!                   **************
!
     &(f,fc,h,fxmat,fxmatpar,
     & volu2d,unsv2d,dt,fxbor,fxborpar,t7,fbor,smh,yasmh,fscexp,
     & nseg,npoin,nptfr,gloseg,sizglo,nbor,limtra,kdir,kddl,optsou,hlin,
     & iopt2,flbortra,surnit,mesh,rain,pluie,train,massou,
     & mass_balance)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    COMPUTES THE TRACER FOR FINITE VOLUME SCHEME.
!+                TO COMPLETE.
!
!history  C-T PHAM (LNHE)
!+        27/02/09
!+        V5P9
!+   JMH : DISTINGUISHES BETWEEN FXBOR AND FXBORTRA
!
!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)
!+        24/02/2012
!+        V6P2
!+   Rain and evaporation added (after initiative by O. Boutron, from
!+   Tour du Valat, and O. Bertrand, Artelia group)
!
!history  J-M HERVOUET (LNHE)
!+        08/06/2015
!+        V7P1
!+   Now mass balance done on request.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| 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
!| FXBORPAR       |-->| FLUXES ON BOUNDARIES (DEFINED ON ALL DOMAIN
!|                |   | AND ASSEMBLED IN PARALLEL)
!| FXMAT          |-->| FLUXES (NON ASSEMBLED IN PARALLEL)
!| FXMATPAR       |-->| FLUXES (ASSEMBLED IN PARALLEL)
!| GLOSEG         |-->| GLOBAL NUMBER OF THE 2 POINTS OF A SEGMENT.
!| H              |-->| WATER DEPTH AT TIME N+1
!|                |   | (VALUE COMPATIBLE WITH CONTINUITY)
!| HLIN           |-->| WATER DEPTH AT TIME N+1
!|                |   | (WITH LINEAR INTERPOLATION IN TIME BETWEEN
!|                |   | HN AND H)
!| 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
!| MASS_BALANCE   |-->| IF YES, ALL TERMS FOR MASS BALANCE
!|                |   | WILL BE COMPUTED
!| MASSOU         |-->| MASS ADDED BY SOURCE TERM
!| MESH           |-->| MESH STRUCTURE
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| 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
!| SIZGLO         |-->| FIRST DIMENSION OF GLOSEG
!| SMH            |-->| SOURCE TERM IN CONTINUITY EQUATION
!| SURNIT         |-->| SURNIT=1/NIT
!| T7             |<->| BIEF_OBJ STRUCTURE FOR A WORK ARRAY
!| TRAIN          |-->| VALUE OF TRACER IN RAIN
!| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
!| YASMH          |-->| IF YES, SMH MUST BE TAKEN INTO ACCOUNT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_tvf => tvf
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NSEG,NPOIN,NPTFR,KDIR,KDDL
      INTEGER, INTENT(IN)             :: SIZGLO,OPTSOU,IOPT2
      INTEGER, INTENT(IN)             :: GLOSEG(sizglo,2)
      INTEGER, INTENT(IN)             :: NBOR(nptfr),LIMTRA(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DT,SURNIT,TRAIN
      DOUBLE PRECISION, INTENT(INOUT) :: FLBORTRA(nptfr)
      DOUBLE PRECISION, INTENT(INOUT) :: F(npoin),MASSOU
      DOUBLE PRECISION, INTENT(IN)    :: FXBOR(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: FC(npoin),H(npoin),HLIN(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: SMH(*),UNSV2D(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: PLUIE(*),VOLU2D(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: FSCEXP(*)
      DOUBLE PRECISION, INTENT(IN)    :: FBOR(nptfr),FXBORPAR(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: FXMAT(nseg),FXMATPAR(nseg)
      LOGICAL, INTENT(IN)             :: YASMH,RAIN,MASS_BALANCE
      TYPE(bief_obj), INTENT(INOUT)   :: T7
      TYPE(bief_mesh), INTENT(INOUT)  :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,N
!
!-----------------------------------------------------------------------
!
      IF(iopt2.EQ.0) THEN
!       CONSERVATIVE ADVECTION FIELD
        DO i = 1,npoin
          f(i) = fc(i)
        ENDDO
      ELSEIF(iopt2.EQ.1) THEN
!       NON CONSERVATIVE ADVECTION FIELD
        DO i = 1,npoin
          f(i) = fc(i)*max(h(i),1.d-8)/max(hlin(i),1.d-8)
        ENDDO
      ELSE
        WRITE(lu,*) 'TVF: UNKNOWN OPTION'
        CALL plante(1)
        stop
      ENDIF
!
!     THE CONTRIBUTION OF FLUXES IS BUILT APART FOR
!     PRELIMINARY PARALLEL ASSEMBLING BEFORE ADDING ON F
!
      DO i = 1,npoin
        t7%R(i) = 0.d0
      ENDDO
      DO i = 1,nseg
        IF(fxmatpar(i).LT.0.d0) THEN
          t7%R(gloseg(i,1)) = t7%R(gloseg(i,1))
     &    - dt/hlin(gloseg(i,1))*unsv2d(gloseg(i,1))
     &    *fxmat(i)*(fc(gloseg(i,2))-fc(gloseg(i,1)))
        ELSEIF(fxmatpar(i).GT.0.d0) THEN
          t7%R(gloseg(i,2)) = t7%R(gloseg(i,2))
     &    + dt/hlin(gloseg(i,2))*unsv2d(gloseg(i,2))
     &    *fxmat(i)*(fc(gloseg(i,1))-fc(gloseg(i,2)))
        ENDIF
      ENDDO
      IF(ncsize.GT.1) CALL parcom(t7,2,mesh)
      DO i = 1,npoin
        f(i) = f(i)+t7%R(i)
      ENDDO
!
!     SOURCE TERMS
!
      IF(yasmh) THEN
        IF(optsou.EQ.1) THEN
          DO i=1,npoin
            IF(smh(i).GT.0.d0) THEN
              f(i)=f(i)+dt/hlin(i)*smh(i)*(fscexp(i)-fc(i))
            ENDIF
          ENDDO
        ELSEIF(optsou.EQ.2) THEN
          DO i=1,npoin
            IF(smh(i).GT.0.d0) THEN
              f(i)=f(i)+dt/hlin(i)*unsv2d(i)*smh(i)*(fscexp(i)-fc(i))
            ENDIF
          ENDDO
        ENDIF
        IF(mass_balance) THEN
          IF(optsou.EQ.1) THEN
            DO i=1,npoin
              IF(smh(i).GT.0.d0) THEN
                massou=massou+dt*volu2d(i)*smh(i)*fscexp(i)
              ELSE
                massou=massou+dt*volu2d(i)*smh(i)*fc(i)
              ENDIF
            ENDDO
          ELSEIF(optsou.EQ.2) THEN
            DO i=1,npoin
              IF(smh(i).GT.0.d0) THEN
                massou=massou+dt*smh(i)*fscexp(i)
              ELSE
                massou=massou+dt*smh(i)*fc(i)
              ENDIF
            ENDDO
          ENDIF
        ENDIF
      ENDIF
!
!     RAIN-EVAPORATION
!
      IF(rain) THEN
        DO i=1,npoin
          IF(pluie(i).GT.0.d0) THEN
!           REAL RAIN, VALUE IN RAIN CONSIDERED...
            f(i)=f(i)+dt/hlin(i)*pluie(i)*(train-fc(i))
          ELSE
!           EVAPORATION, VALUE IN RAIN NOT CONSIDERED...
            f(i)=f(i)+dt/hlin(i)*pluie(i)*(     -fc(i))
          ENDIF
        ENDDO
      ENDIF
!
!     ON THE DIRICHLET BOUNDARIES, FLUX TERMS TAKEN INTO ACCOUNT
!     ON OTHERS, FBOR IS TAKEN AS FN, SO NO CONTRIBUTION
!
      DO i=1,nptfr
        IF(limtra(i).EQ.kdir) THEN
          n=nbor(i)
          f(n)=f(n)-dt/hlin(n)*unsv2d(n)*fxborpar(n)*(fbor(i)-fc(n))
        ENDIF
      ENDDO
      IF(mass_balance) THEN
        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)*fc(n)*surnit
          ENDIF
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END