cflp12.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE cflp12
!                   *****************
!
     &(u,v,x,y,ikle,nelem,nelmax,w1)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE COURANT NUMBER AT EACH POINT OF THE MESH
!+                AND FOR EACH TIMESTEP.
!+
!+            THE STABILITY CRITERION OF THE DISTRIBUTIVE SCHEME N
!+                IS HERE USED TO EVALUATE THE COURANT NUMBER.
!
!history  C MOULIN   (LNH)
!+
!+
!+
!
!history  JMH
!+        29/12/05
!+        V5P6
!+   MODIFICATIONS
!
!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  S.E.BOURBAN (HRW)
!+        21/03/2017
!+        V7P3
!+   Replacement of the DATA declarations by the PARAMETER associates
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IKLE           |-->| CONNECTIVITY TABLE
!| NELEM          |-->| NUMBER OF ELEMENTS IN THE MESH
!| NELMAX         |-->| FIRST DIMENSION OF IKLE, MAXIMUM NUMBER OF ELEMENTS
!|                |   | IN THE MESH
!| U              |-->| VELOCITY ALONG X.
!| V              |-->| VELOCITY ALONG Y.
!| W1             |-->| RESULT IN NON ASSEMBLED FORM
!| X              |-->| ABSCISSAE OF POINTS GIVEN PER ELEMENT
!| Y              |-->| ORDINATES OF POINTS GIVEN PER ELEMENT
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)  :: NELEM,NELMAX
      DOUBLE PRECISION, INTENT(IN)  :: U(*),V(*)
      DOUBLE PRECISION, INTENT(IN)  :: X(nelmax*3),Y(nelmax*3)
      INTEGER         , INTENT(IN)  :: IKLE(nelmax*4)
      DOUBLE PRECISION, INTENT(OUT) :: W1(nelmax*4)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,IT,IAD1,IAD2,IAD3,IG1,IG2,IG3
!
      DOUBLE PRECISION USUR2,VSUR2
      DOUBLE PRECISION SUR6,K1,K2,K3,L12,L13,L21,L23,L31,L32
      DOUBLE PRECISION X1,X2,X3,Y1,Y2,Y3,TIERS
!
      INTRINSIC max,min
!
!-----------------------------------------------------------------------
!
!     FOR A QUASI-BUBBLE TRIANGLE : NUMBERS OF THE VERTICES OF THE
!     SUB-TRIANGLES IN THE INITIAL TRIANGLE
!     IL(NUMBER OF THE SUB-TRIANGLE,LOCAL NUMBER IN THE SUB-TRIANGLE)
!
      INTEGER :: IL(3,3)
      parameter( il = reshape( (/
     &          1,2,3,2,3,1,4,4,4 /), shape=(/ 3,3 /) ) )
!
!-----------------------------------------------------------------------
!
      tiers= 1.d0 / 3.d0
      sur6 = 1.d0 / 6.d0
!
!     INITIALISES W
!
      DO ielem = 1 , 4*nelmax
        w1(ielem) = 0.d0
      ENDDO ! IELEM
!
!     USING THE PSI SCHEME,
!     LOOP ON THE 3 SUB-TRIANGLES AND PRE-ASSEMBLY
!
      DO it=1,3
      DO ielem = 1 , nelem
!
!       ADDRESSES IN AN ARRAY (NELMAX,*)
        iad1= ielem + (il(it,1)-1)*nelmax
        iad2= ielem + (il(it,2)-1)*nelmax
        iad3= ielem + (il(it,3)-1)*nelmax
!       GLOBAL NUMBERS IN THE INITIAL TRIANGLE
        ig1 = ikle(iad1)
        ig2 = ikle(iad2)
        ig3 = ikle(iad3)
!       COORDINATES OF THE SUB-TRIANGLE VERTICES
        x1 = x(iad1)
        x2 = x(iad2) - x1
        y1 = y(iad1)
        y2 = y(iad2) - y1
!       POINT 3 IS ALWAYS AT THE CENTRE OF THE INITIAL TRIANGLE
        x3=tiers*(x(ielem)+x(ielem+nelmax)+x(ielem+2*nelmax))-x1
        y3=tiers*(y(ielem)+y(ielem+nelmax)+y(ielem+2*nelmax))-y1
!
        usur2 = (u(ig1)+u(ig2)+u(ig3))*sur6
        vsur2 = (v(ig1)+v(ig2)+v(ig3))*sur6
!
        k1 = usur2 * (y2-y3) - vsur2 * (x2-x3)
        k2 = usur2 * (y3   ) - vsur2 * (x3   )
        k3 = usur2 * (  -y2) - vsur2 * (  -x2)
!
        l12 = max(  min(k1,-k2) , 0.d0 )
        l13 = max(  min(k1,-k3) , 0.d0 )
        l21 = max(  min(k2,-k1) , 0.d0 )
        l23 = max(  min(k2,-k3) , 0.d0 )
        l31 = max(  min(k3,-k1) , 0.d0 )
        l32 = max(  min(k3,-k2) , 0.d0 )
!
        w1(iad1) = w1(iad1) + l12 + l13
        w1(iad2) = w1(iad2) + l21 + l23
        w1(iad3) = w1(iad3) + l31 + l32
!
      ENDDO ! IELEM
      ENDDO ! IT
!
!-----------------------------------------------------------------------
!
      RETURN
      END