pentco.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE pentco
!                       *****************
     &(ii)
!
!***********************************************************************
! ARTEMIS   V7P3                                     Aug 2017
!***********************************************************************
!
!
!brief FUNCTION:   CALCULATE SECOND ORDER BOTTOM EFFECTS (GRADIENT&CURVATURE)
!+                  FOR EXTENDED MILD-SLOPE EQUATION
!
!code
!+        OUTPUT :
!+        1 + F = 1 + E1(KH)*grad(H)**2 + E2(KH)/K0*LAPLACIEN(H)
!+
!+      WE CAN CHOOSE TO ONTEGRATE ONLY GRADIANT EFFECTS              : II=1
!+      WE CAN CHOOSE TO ONTEGRATE ONLY CURVATURE EFFECTS             : II=2
!+      WE CAN CHOOSE TO INTEGRATE BOTH GRADIANT & CURVATURE EFFECTS  : II=3
!+      DEFAULT VALUE  : IPENTCO=0, ARTEMIS SOLVE CLASSICAL MILD-SLOPE EQUATION
!+
!+
!+     EXPRESSIONS FOR E1 and  E2 USED HERE ARE (Chamberlain & Porter 1995) :
!+
!+     (given X = 2 KH)
!+
!+          ( X**4 + 4 X**3 SH(X) - 9 SH(X)SH(2X) + 3 X (X+2SH(X))*(CH(X)**2-2CH(X)+3) )
!+ E1(KH) = -----------------------------------------------------------------------------
!+                                  3 ( X+SH(X) )**4
!+
!+
!+
!+
!+                         (TH(X)-X)* CH(X)
!+    E2(KH)/K0 = 2 H * ------------------------
!+                         X ( SH(X) + X )**2
!+
!+      K0 IS THE WAVE NUMBER FOR INFINITE DEPTH : K0 = K TH(KH)
!+      H IS THE WATER DEPTH
!+                                           ------
!+
!+  USING SUBROUTINE BERKHO NOTATIONS,
!+  AFTER VARIATIONAL FORMULATION :
!+           /
!+ AM1 =    / C*CG * GRAD(PSII)*GRAD(PSIJ) DS
!+         /S
!+
!+           /
!+       -  / OMEGA**2 * CG/C * (1+F) * PSII*PSIJ  DS
!+         /S
!+
!+           /
!+       -  /  BPHIRB * PSII*PSIJ  DB
!+         /B
!+
!+
!+ THE SECOND MEMEBER (DIFFUSION) IS MODIFIED
!
!history C.PEYRARD & E.RAZAFINDRAKOTO
!+        31/05/11
!+        V6P1
!+
!
!history  N.DURAND (HRW)
!+        August 2017
!+        V7P3
!+   Removed unnecessary references to PI, DEGRAD and RADDEG
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!|    II          |-->|  OPTION FOR GRADIENT AND CURVATURE EFFECTS
!|    T2          |---|  WORK TABLE
!|    T4          |---|  WORK TABLE
!|    T5          |---|  WORK TABLE
!|    T6          |---|  WORK TABLE
!|    T7          |---|  WORK TABLE
!|    T9          |---|  WORK TABLE
!|    T8          |---|  WORK TABLE
!|    T11         |---|  WORK TABLE
!|    T12         |---|  WORK TABLE
!|    T3          |<--|  OUTPUT WORK TABLE WITH CORRECTION TERMS FOR
!|                |   |  GRADIENT AND CURVATURE EFFECTS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! APPELE PAR                   :  BERKHO
! SOUS-PROGRAMMES APPELES      :  FCTE1 et FCTE2
! TABLEAUX DE TRAVAIL UTILISES :  T2 T3 T4 T5 T6 T7 T9 T8 T11 T12
!
!-----------------------------------------------------------------------
!
      USE bief
      USE interface_artemis,ONLY: fcte1,fcte2
      USE declarations_telemac
      USE declarations_artemis
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: II
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I
!
!-----------------------------------------------------------------------
!
!
! MASS MATRIX
!
      CALL vector(t2 , '=' , 'MASBAS          ' , ielm ,
     &            1.d0 , c , c , c , c , c , c ,
     &            mesh , msk  , maskel )
!
      IF(ii.EQ.1.OR.ii.EQ.3) THEN
! --------------
! GRADIENT EFFECTS
! --------------
!
!--->  E1*GRAD(H)**2 ---> IN T4
!
!  DX
!
      CALL vector(t7 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , h , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
      CALL os('X=YZ    ', x=t4, y=t7, z=t7)
!
!  DY
!
      CALL vector(t7 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , h , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
      CALL os('X=YZ    ' , x=t5,y=t7,z=t7)
!
! GRAD(H)**2
!
      CALL os( 'X=X+Y   ', x=t4,y=t5)
!
! SQUARE MASS
!
      CALL os( 'X=YZ    ' , x=t8,y=t2,z=t2)
!
! COEFF GRAD(H)**2
!
      CALL os( 'X=Y/Z   ' , x=t4,y=t4,z=t8)
!
!---> FUNCTION E1
!
      DO i=1,npoin
        t11%R(i)=fcte1(  k%R(i)*h%R(i) )
      END DO
!
!--->  E1*GRAD(H)**2 ---> IN T4
!
      CALL os( 'X=YZ    ' , x=t4,y=t4,z=t11)
! END OF GRADIENT EFFECTS
      ENDIF
!
      IF(ii.EQ.2.OR.ii.EQ.3) THEN
! ------------------
! CURVATURE EFFECTS
! ------------------
!---> E2/K0*LAPLACIAN(H)   ---> IN T9
!  DX
      CALL vector(t7 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , h , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
      CALL os( 'X=Y/Z    ' , x=t7,y=t7,z=t2)

      CALL vector(t5 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , t7 , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
!
! COEFF DX
!
      CALL os('X=Y/Z    ' , x=t5,y=t5,z=t2)
!
!  DY
!
      CALL vector(t7 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , h , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
      CALL os('X=Y/Z    ' , x=t7,y=t7,z=t2)
!
      CALL vector(t6 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , t7 , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
!
! COEFF DY
!
      CALL os('X=Y/Z   ', x=t6,y=t6,z=t2)
!
      CALL os('X=Y+Z   ', x=t9,y=t5,z=t6)
!
!---> FUNCTION E2 * 2 H
!
      DO i=1,npoin
        t12%R(i)=2.*h%R(i) * fcte2(k%R(i)*h%R(i))
      END DO
!
!---> E2/K0*LAPLACIAN(H)
!
      CALL os( 'X=YZ    ' , x=t9,y=t9,z=t12)
!
! END OF CURVATURE EFFECTS
      ENDIF
!
! SUM OF GRADIENT AND CURVTURE EFFECTS, DEPENDING OF OPTION "IPENTCO"
!
      IF(ii.EQ.1)  THEN
!       F= E1*GRAD(H)**2
        CALL os( 'X=Y      ' , x=t3,y=t4)
      ENDIF
!
      IF(ii.EQ.2)  THEN
!       F= E2/K0*LAPLACIAN(H)
        CALL os( 'X=Y      ' , x=t3,y=t9)
      ENDIF
!
      IF(ii.EQ.3)  THEN
!       F =  E1*grad(H)**2 + E2/K0*LAPLACIEN(H)
        CALL os( 'X=Y+Z    ' , x=t3,y=t4,z=t9)
      ENDIF
!
! ADD 1.,  T3 = 1 + F
!
      CALL os( 'X=X+C   ', x=t3,c=1.d0)
!
!-----------------------------------------------------------------------
!
      RETURN
      END