flux_ef_vf_3d.f

Go to the documentation of this file.

!                   ************************
                    SUBROUTINE flux_ef_vf_3d
!                   ************************
!
     &(flow,w2d,w3d,nseg2d,nelem2,nelmax2,mesh2d,iniflo,
     & iopt,sens,ielm3,nplan,ikle,nelmax,knolg)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    The 3D Element By Element fluxes leaving points are grouped
!+        on every plane (one contribution from below, one from above).
!+        Then every plane is considered like a 2D mesh, and the 2D
!+        EBE fluxes are transformed into horizontal fluxes between
!+        points (with a call to flux_ef_vf). This gives all the 3D
!+        horizontal fluxes between points.
!+
!+
!warning  For element 51 the result is summed on the vertical, and is
!+        a 2D result (see correction_depth_3d where flux_ef_vf_3d is
!+        called, otherwise this routine should not be called for
!+        tetrahedra).
!+
!
!history
!+        15/05/2009
!+        V6P0
!+   INSPIRED FROM LEO POSTMA (DELTARES)
!
!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)
!+        26/08/2011
!+        V6P2
!+   Adaptation to element 51 (prisms cut into tetrahedrons)
!
!history  J.-M. HERVOUET (LNHE)
!+        24/08/2012
!+        V6P2
!+   Computation of fluxes and assembly swapped in intermediate planes.
!
!history  J.-M. HERVOUET (LNHE)
!+        11/09/2017
!+        V7P3
!+   Argument NELMAX2 added for dimensioning arrays.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FLOW           |<--| FLUX
!| IELM3          |-->| DISCRETISATION IN 3D
!|                |   | (41: PRISMS 51:PRISMS CUT INTO TETRAHEDRONS)
!| IKLE           |-->| CONNECTIVITY TABLE
!| INIFLO         |-->| IF(YES) FLOW WILL BE INITIALISED AT 0.
!| IOPT           |-->| CHOICE OF THE CONSTANT IN FLUX_EF_VF
!| KNOLG          |-->| GIVES THE ORIGINAL GLOBAL NUMBER OF POINTS
!|                |   | IN SCALAR MODE (SIZE NPOIN3 BUT FILLED ONLY
!|                |   | UP TO NPOIN2)
!| MESH2D         |-->| 2D MESH
!| NELEM2         |-->| NUMBER OF ELEMENTS IN 2D
!| NELMAX2        |-->| MAXIMUM NUMBER OF ELEMENTS IN 2D
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS (SEE IKLE)
!| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH
!| NSEG2D         |-->| NUMBER OF SEGMENTS IN 2D
!| SENS           |-->| IF 1: HORIZONTAL FLUXES FROM BOTTOM TO TOP
!|                |   | IF 2: HORIZONTAL FLUXES FROM TOP TO BOTTOM
!| W2D            |<--| NON ASSEMBLED FLUXES LEAVING POINTS,PER TRIANGLE
!| W3D            |-->| NON ASSEMBLED FLUXES LEAVING POINTS,PER PRISM
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_flux_ef_vf_3d => flux_ef_vf_3d
      USE declarations_telemac, ONLY : tetra
!
      USE declarations_special
      IMPLICIT NONE
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      INTEGER, INTENT(IN)             :: NSEG2D,NELEM2
!                                             *=NSEG2D*NPLAN+NPOIN2*NETAGE
      INTEGER, INTENT(IN)             :: IOPT,SENS,IELM3,NPLAN,NELMAX
      INTEGER, INTENT(IN)             :: NELMAX2
!                                                    6 IF IELM3=41
!                                                    4 IF IELM3=51
      INTEGER, INTENT(IN)             :: IKLE(nelmax,*),KNOLG(*)
      DOUBLE PRECISION, INTENT(INOUT) :: FLOW(*)
!                                                   6 IF IELM3=41
!                                                   4 IF IELM3=51
      DOUBLE PRECISION, INTENT(IN)    :: W3D(nelmax,*)
      DOUBLE PRECISION, INTENT(INOUT) :: W2D(nelmax2,3)
      LOGICAL, INTENT(IN)             :: INIFLO
      TYPE(bief_mesh), INTENT(INOUT)  :: MESH2D
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      INTEGER IPLAN,N1,N2,N3,N4,IELEM,I1,I2,I3,S1,S2,S3,IT,IELEM3D,K,L
!
!     TETRA : WILL GIVE THE LOCAL NUMBERS OF POINTS IN THE PRISM
!             THE 0 CORRESPOND TO SITUATIONS
!             THAT NEVER HAPPEN (TETRA(1,1,1,... OR TETRA(2,2,2,...)
!     INTEGER TETRA(2,2,2,3,4)
!     DATA TETRA / 0,1,1,1,1,1,1,0,0,4,4,4,4,4,4,0,0,6,4,5,5,4,6,0,
!    &             0,2,2,2,2,2,2,0,0,6,6,6,6,6,6,0,0,3,1,2,2,1,3,0,
!    &             0,3,3,3,3,3,3,0,0,5,5,5,5,5,5,0,0,2,3,4,1,6,5,0,
!    &             0,4,5,4,6,6,5,0,0,2,3,3,1,2,1,0,0,4,5,3,6,2,1,0 /
!
!-----------------------------------------------------------------------
!
      IF(sens.NE.1.AND.sens.NE.2) THEN
        WRITE(lu,*) 'FLUX_EF_VF_3D: UNKNOWN SENS = ',sens
        CALL plante(1)
        stop
      ENDIF
!
      IF(sens.EQ.2.AND..NOT.iniflo) THEN
        WRITE(lu,*) 'FLUX_EF_VF_3D: SENS = 2 AND INIFLO = .FALSE.'
        WRITE(lu,*) '               INCOMPATIBLE OPTIONS'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(ielm3.EQ.41) THEN
!
!     CASE OF PRISMS
!
!     ADDS FLUXES ON HORIZONTAL SEGMENTS (FOR SEGMENT NUMBERING IN 3D
!                                         SEE STOSEG41 IN BIEF)
!
!     FOR INTERMEDIATE PLANES
!     HERE WE HAVE THE CHOICE OF COMPUTING THE FLUXES BEFORE OR AFTER
!     ASSEMBLING ON VERTICAL. UP TO VERSION 6.1, ASSEMBLY WAS DONE FIRST
!     FROM VERSION 6.2 ON WE COMPUTE THE FLUXES FIRST AND ASSEMBLE AFTER
!     THIS IS MORE COMPATIBLE WITH DISTRIBUTIVE SCHEMES THAT COMPUTE
!     THEIR FLUXES AT ELEMENT LEVEL. THEN THE FLUX LIMITATION DONE
!     IN NA_FLUX3D_LIM WILL BE COMPATIBLE.
!
      DO iplan=1,nplan
!
        IF(sens.EQ.1) THEN
          n1=(iplan-1    )*nseg2d+1
          n2= iplan       *nseg2d
        ELSE
          n1=(  nplan-iplan)*nseg2d+1
          n2=(1+nplan-iplan)*nseg2d
        ENDIF
!
!       POINTS 1, 2 AND 3 OF UPPER LEVEL
        IF(iplan.EQ.1) THEN
!         FIRST PLANE: ONLY POINTS 1, 2 AND 3 OF UPPER LEVEL
          DO ielem=1,nelem2
            w2d(ielem,1)=w3d(ielem,1)
            w2d(ielem,2)=w3d(ielem,2)
            w2d(ielem,3)=w3d(ielem,3)
          ENDDO
          CALL flux_ef_vf(flow(n1:n2),w2d,nseg2d,nelem2,mesh2d%NELMAX,
     &                    mesh2d%ELTSEG%I,mesh2d%ORISEG%I,
     &                    mesh2d%IKLE%I,iniflo,iopt)
        ELSEIF(iplan.EQ.nplan) THEN
!         LAST PLANE: ONLY POINTS 4, 5 AND 6 OF LOWER LEVEL
          n3=nelem2*(iplan-2)
          DO ielem=1,nelem2
            w2d(ielem,1)=w3d(n3+ielem,4)
            w2d(ielem,2)=w3d(n3+ielem,5)
            w2d(ielem,3)=w3d(n3+ielem,6)
          ENDDO
          CALL flux_ef_vf(flow(n1:n2),w2d,nseg2d,nelem2,mesh2d%NELMAX,
     &                    mesh2d%ELTSEG%I,mesh2d%ORISEG%I,
     &                    mesh2d%IKLE%I,iniflo,iopt)
        ELSE
!       INTERMEDIATE PLANE
!         POINTS 4, 5 AND 6 OF LOWER LEVEL
          n3=nelem2*(iplan-2)
          DO ielem=1,nelem2
            w2d(ielem,1)=w3d(n3+ielem,4)
            w2d(ielem,2)=w3d(n3+ielem,5)
            w2d(ielem,3)=w3d(n3+ielem,6)
          ENDDO
          CALL flux_ef_vf(flow(n1:n2),w2d,nseg2d,nelem2,mesh2d%NELMAX,
     &                    mesh2d%ELTSEG%I,mesh2d%ORISEG%I,
     &                    mesh2d%IKLE%I,iniflo,iopt)
!         POINTS 1, 2, 3 OF UPPER LEVEL
          n4=n3+nelem2
          DO ielem=1,nelem2
            w2d(ielem,1)=w3d(n4+ielem,1)
            w2d(ielem,2)=w3d(n4+ielem,2)
            w2d(ielem,3)=w3d(n4+ielem,3)
          ENDDO
          CALL flux_ef_vf(flow(n1:n2),w2d,nseg2d,nelem2,mesh2d%NELMAX,
     &                    mesh2d%ELTSEG%I,mesh2d%ORISEG%I,
     &                    mesh2d%IKLE%I,.false.,iopt)
!                                       !!!!!!!
        ENDIF
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      ELSEIF(ielm3.EQ.51) THEN
!
!     CASE OF PRISMS CUT INTO TETRAHEDRONS
!
!     ADDS FLUXES ON HORIZONTAL SEGMENTS (FOR SEGMENT NUMBERING IN 3D
!                                         SEE STOSEG51 IN BIEF)
!
!
!     INITIALISING
!
      DO ielem=1,nelem2
        w2d(ielem,1)=0.d0
        w2d(ielem,2)=0.d0
        w2d(ielem,3)=0.d0
      ENDDO
!
!     LOOP ON PLANES
!
      DO iplan=1,nplan
!
!       LOOP ON TRIANGLES
!
        DO ielem=1,nelem2
!
!         HERE LOWER LEVEL OF ELEMENTS ALWAYS TAKEN
!         TO FIND THE WAY THE PRISM HAS BEEN CUT BY LOOKING
!         AT GLOBAL NUMBERS OF POINTS
!         THIS PART IS THUS COMMON TO ALL PLANES
!         IKLE 3D IS TAKEN, COULD BE IKLE 2D AS WELL
!         SO NO DEPENDENCE ON LEVEL IPLAN
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
!         IN PARALLELISM BACK TO ORIGINAL GLOBAL NUMBERS
          IF(ncsize.GT.1) THEN
            i1=knolg(i1)
            i2=knolg(i2)
            i3=knolg(i3)
          ENDIF
!         THIS IS DONE LIKE IN CPIKLE3 TO USE ARRAY TETRA
          IF(i1.GT.i2) THEN
            s1=1
          ELSE
            s1=2
          ENDIF
          IF(i2.GT.i3) THEN
            s2=1
          ELSE
            s2=2
          ENDIF
          IF(i3.GT.i1) THEN
            s3=1
          ELSE
            s3=2
          ENDIF
!
!         NOW TAKING CONTRIBUTIONS OF TETRAHEDRON K= 1, 2 AND 3
!         OF LOWER AND UPPER LEVEL
!
!         LOWER LEVEL : POINTS 4,5,6 OF ORIGINAL PRISM CONTRIBUTE
!
          IF(iplan.GT.1) THEN
!           LOOP ON 3 TETRAHEDRONS
            DO k=1,3
              ielem3d=3*(iplan-2)*nelem2+ielem+(k-1)*nelem2
!             POINTS 1 TO 4
              DO l=1,4
!               ORIGINAL NUMBER IN THE PRISM
                it=tetra(s1,s2,s3,k,l)
                IF(it.GE.4) THEN
                  w2d(ielem,it-3)=w2d(ielem,it-3)+w3d(ielem3d,l)
                ENDIF
              ENDDO
            ENDDO
          ENDIF
!
!         UPPER LEVEL : POINTS 1,2,3 OF ORIGINAL PRISM CONTRIBUTE
!
          IF(iplan.LT.nplan) THEN
!           LOOP ON 3 TETRAHEDRONS
            DO k=1,3
              ielem3d=3*(iplan-1)*nelem2+ielem+(k-1)*nelem2
!             POINTS 1 TO 4
              DO l=1,4
!               ORIGINAL NUMBER IN THE PRISM
                it=tetra(s1,s2,s3,k,l)
                IF(it.LE.3) THEN
                  w2d(ielem,it)=w2d(ielem,it)+w3d(ielem3d,l)
                ENDIF
              ENDDO
            ENDDO
          ENDIF
!
        ENDDO
!
        IF(sens.EQ.1) THEN
          n1=(iplan-1    )*nseg2d+1
          n2= iplan       *nseg2d
        ELSE
          n1=(  nplan-iplan)*nseg2d+1
          n2=(1+nplan-iplan)*nseg2d
        ENDIF
!
      ENDDO
!
!     HERE 2D FLUXES !!!!!!!!!!!!!
!
      CALL flux_ef_vf(flow,w2d,nseg2d,nelem2,mesh2d%NELMAX,
     &                mesh2d%ELTSEG%I,mesh2d%ORISEG%I,
     &                mesh2d%IKLE%I,iniflo,iopt)
!
!-----------------------------------------------------------------------
!
      ELSE
        WRITE(lu,*) 'UNEXPECTED ELEMENT TYPE IN FLUX_EF_VF_3D'
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END