mt15pp.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE mt15pp
!                   *****************
!
     &(t,xm,xmul,f,zpt,surfac,ikle,nelem,nelmax)
!
!***********************************************************************
! BIEF   V7P0                                  03/06/2014
!***********************************************************************
!
!brief    Builds a matrix corresponding to the advection with a settling
!+        velocity. This vertical velocity WC is stored in F and
!+        is positive when going downwards.
!
!warning  Element linear prism : 41
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        03/06/2014
!+        V7P0
!+   First version.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| F              |-->| FUNCTION USED IN THE FORMULA
!| IKLE           |-->| CONNECTIVITY TABLE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| SURFAC         |-->| AREA OF 2D ELEMENTS
!| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
!| XM             |<->| OFF-DIAGONAL TERMS
!| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
!| ZPT            |-->| Z COORDINATES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_mt15pp => mt15pp
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX
      INTEGER, INTENT(IN) :: IKLE(nelmax,6)
!
      DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(nelmax,30)
!
      DOUBLE PRECISION, INTENT(IN) :: XMUL
      DOUBLE PRECISION, INTENT(IN) :: F(*),ZPT(*)
!
      DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
!
      INTEGER IELEM,I1,I2,I3,I4,I5,I6
      DOUBLE PRECISION XSUR3
!
!-----------------------------------------------------------------------
!
      xsur3=xmul/3.d0
!
!=======================================================================
!
      DO ielem=1,nelem
!
        i1=ikle(ielem,1)
        i2=ikle(ielem,2)
        i3=ikle(ielem,3)
        i4=ikle(ielem,4)
        i5=ikle(ielem,5)
        i6=ikle(ielem,6)
!
        t(ielem,1)=0.d0
        t(ielem,2)=0.d0
        t(ielem,3)=0.d0
!
        xm(ielem, 1)=0.d0
        xm(ielem, 2)=0.d0
        xm(ielem, 4)=0.d0
        xm(ielem, 5)=0.d0
        xm(ielem, 6)=0.d0
        xm(ielem, 7)=0.d0
        xm(ielem, 9)=0.d0
        xm(ielem,10)=0.d0
        xm(ielem,11)=0.d0
        xm(ielem,13)=0.d0
        xm(ielem,14)=0.d0
        xm(ielem,15)=0.d0
!
        xm(ielem,16)=0.d0
        xm(ielem,17)=0.d0
        xm(ielem,18)=0.d0
        xm(ielem,19)=0.d0
        xm(ielem,20)=0.d0
        xm(ielem,21)=0.d0
        xm(ielem,22)=0.d0
        xm(ielem,23)=0.d0
        xm(ielem,24)=0.d0
        xm(ielem,25)=0.d0
        xm(ielem,26)=0.d0
        xm(ielem,27)=0.d0
        xm(ielem,28)=0.d0
        xm(ielem,29)=0.d0
        xm(ielem,30)=0.d0
!
!-----------------------------------------------------------------------
!
!     EXTRA-DIAGONAL TERMS: BOTTOM POINTS RECEIVE FROM UPPER POINT
!                           A FLUX EQUAL TO WC*SURFAC/3 MULTIPLIED
!                           BY CONCENTRATION OF UPPER POINT
!     DIAGONAL TERMS : TOP POINTS GIVE TO LOWER POINTS
!                      A FLUX EQUAL TO WC*SURFAC/3 MULTIPLIED
!                      BY ITS CONCENTRATION.
!
!-----------------------------------------------------------------------
!
        IF(zpt(i4)-zpt(i1).GT.1.d-4) THEN
!         TERM 4-4
          t(ielem,4)=   f(i4)*surfac(ielem)*xsur3
!         TERM 1-4
          xm(ielem, 3)=-f(i4)*surfac(ielem)*xsur3
        ENDIF
        IF(zpt(i5)-zpt(i2).GT.1.d-4) THEN
!         TERM 5-5
          t(ielem,5)=   f(i5)*surfac(ielem)*xsur3
!         TERM 2-5
          xm(ielem, 8)=-f(i5)*surfac(ielem)*xsur3
        ENDIF
        IF(zpt(i6)-zpt(i3).GT.1.d-4) THEN
!         TERM 6-6
          t(ielem,6)=   f(i6)*surfac(ielem)*xsur3
!         TERM 3-6
          xm(ielem,12)=-f(i6)*surfac(ielem)*xsur3
        ENDIF
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END