qporos.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE qporos
!                   *****************

     &( tstot , tsder , f , cg,  lt,xk,
     &  nf    , ndire  , npoin2   , amorp )

!***********************************************************************
! TOMAWAC   V7P3
!***********************************************************************
!
!brief    Takes into account the friction due to porous media
!
!history  THIERRY FOUQUET (EDF - LNHE)
!+        27/03/2018
!+        V7P3
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AMORP          |<->| DUMPING DUE TO POROUS MEDIA
!| CG             |<--| DISCRETIZED GROUP VELOCITY
!| F              |-->| DIRECTIONAL SPECTRUM
!| LT             |-->| TIME STEP
!| NF             |-->| NUMBER OF FREQUENCIES
!| NDIRE          |-->| NUMBER OF DIRECTIONS
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| TSDER          |<->| DERIVED PART OF THE SOURCE TERM CONTRIBUTION
!| TSTOT          |<->| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
!| XK             |<--| WAVE NUMBER
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!

      USE declarations_tomawac, ONLY: deupi, gravit, freq, depth, x, y
!
      USE interface_tomawac, ex_poros => qporos
      USE declarations_special
      IMPLICIT NONE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      INTEGER, INTENT(IN)             :: NF,NDIRE,NPOIN2,LT
      DOUBLE PRECISION, INTENT(INOUT) :: AMORP(npoin2,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: TSTOT(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: TSDER(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: CG(npoin2,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: XK(npoin2,nf)
      DOUBLE PRECISION, INTENT(IN)    :: F(npoin2,ndire,nf)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      INTEGER N
! THe complexe dispersive function and its derivative
      COMPLEX(KIND=R8) FONC, DERFONC, SECONDFONC, DERTROIS
      EXTERNAL fonc, derfonc, secondfonc, dertrois
! ALPHA MIGHT BE READ IN A FILE OR SET BY POLYNOME
      DOUBLE PRECISION GAMMA, ALPHA, INIT
      COMPLEX(KIND=R8) X0,X1,X2,X3, PHI, PHI1
      COMPLEX(KIND=R8) LPRIM,INVLPRIM,LSECOND,LTIERCE
      COMPLEX(KIND=R8) TEMP
      COMPLEX(KIND=R8) FPRIM,FSECOND
      DOUBLE PRECISION TAU,COEFF,CA,S,HW
      INTEGER IP,JF,JP
      IF (lt.EQ.1) THEN
!     TAUX DE POROSITE
        tau=0.8d0
!     COEFFICIENT DE MASSE AJOUTEE
        ca=2.5d0
!     COEFFICIENT D'AMORTISSEMENT
        coeff=1.d0
!     AUTRE COEFF POROSITE
        s=1+ca*(1.d0-tau)/tau
        temp=cmplx(s,1.d0, kind=r8)
!     COEFFICIENT DE POROSITE
        phi=tau*temp/(s**2+coeff**2)
        alpha=0.8d0
        phi1=cmplx(1.d0,0.d0,kind=r8)
        DO ip=1,npoin2
          IF(     x(ip).GE.40.d0.AND.x(ip).LE.50.d0
     &       .AND.y(ip).GE.39.d0 .AND.y(ip).LE.61.d0)
     &          THEN
            alpha=0.8d0
            hw=(1-alpha)*depth(ip)
            DO jf=1,nf
              gamma=(deupi*freq(jf))**2*depth(ip)/gravit
              init=1.d0
              DO n=1,10
                init=gamma/tanh(init)
              ENDDO
              x0=cmplx(init,0.d0, kind=r8)
!      CALCULATION OF COMPLEXE WAVE NUMBER ACCORDING TO Liou 2006
              lprim=derfonc(x0,gamma,phi1,alpha)
              invlprim=1.d0/lprim
              fprim=derfonc(x0,gamma,phi,alpha)
              fsecond=secondfonc(x0,gamma,phi,alpha)
              lsecond=secondfonc(x0,gamma,phi1,alpha)
              ltierce=dertrois(x0, gamma, phi1,alpha)
              x1=-fonc(x0,gamma,phi,alpha)*invlprim
              x2=-invlprim*(0.5d0*lsecond*x1**2+(fprim-lprim)*x1)
              x3=-invlprim*(lsecond*x1*x2+(fprim-lprim)*x2
     &           +0.5d0*(fsecond-lsecond)*x1**2+1.d0/6.d0*x1**3*ltierce)
!     CALCULATION OF CELERITY FOR THE REAL WATER DEPTH
              gamma=(deupi*freq(jf))**2*hw/gravit
              init=deupi*freq(jf)*hw/sqrt(gravit*hw)
              DO n=1,50
                init=gamma/tanh(init)
              ENDDO
              cg(ip,jf)=0.5*(1+2*init/sinh(2*init))
     &             *deupi*freq(jf)*hw/init
              xk(ip,jf)=init/depth(ip)

              amorp(ip,jf)=2.d0*cg(ip,jf)*aimag(x0+x1+x2+x3)
            ENDDO
          ELSE
            DO jf=1,nf
              amorp(ip,jf)=0.d0
            ENDDO
          ENDIF
        ENDDO
      ENDIF                     !LT.LT.1
      DO ip=1,npoin2
        DO jf=1,nf
          DO jp=1,ndire
            tstot(ip,jp,jf) = tstot(ip,jp,jf)+amorp(ip,jf)*f(ip,jp,jf)
            tsder(ip,jp,jf) = tsder(ip,jp,jf)+amorp(ip,jf)
          ENDDO
        ENDDO
      ENDDO
      END

!-----------------------------------------------------------------------
      FUNCTION fonc(X, GAMMA, PHI, ALPHA)
!
      IMPLICIT NONE
      COMPLEX(KIND(1.D0)) FONC
      DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
      COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
      COMPLEX(KIND(1.D0)) TANHCPLX
      EXTERNAL tanhcplx

      DOUBLE PRECISION BETA
      beta=1.d0-alpha
      fonc=gamma-x*tanhcplx(beta*x)
     &          -phi*tanhcplx(alpha*x)*(x-gamma*tanhcplx(beta*x))
      RETURN
      END
!-----------------------------------------------------------------------
      FUNCTION derfonc(X, GAMMA, PHI, ALPHA)
!
      IMPLICIT NONE
      COMPLEX(KIND(1.D0)) DERFONC
      DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
      COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
      COMPLEX(KIND(1.D0)) TANHCPLX,COSHCPLX
      EXTERNAL tanhcplx,coshcplx

      DOUBLE PRECISION BETA
      beta=1.d0-alpha
      derfonc=-phi*tanhcplx(alpha*x)
     &            *(1.d0-gamma*beta/coshcplx(beta*x)**2)
     &        -phi*alpha*(x-gamma*tanhcplx(beta*x))/coshcplx(alpha*x)**2
     &     -beta*x/coshcplx(beta*x)**2-tanhcplx(beta*x)
      RETURN
      END
!-----------------------------------------------------------------------
      FUNCTION secondfonc(X, GAMMA, PHI,  ALPHA)
!
      IMPLICIT NONE
      COMPLEX(KIND(1.D0)) SECONDFONC
      DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
      COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
      COMPLEX(KIND(1.D0)) TANHCPLX,COSHCPLX
      EXTERNAL tanhcplx,coshcplx

      DOUBLE PRECISION BETA
      beta=1-alpha

      secondfonc=
     & -2.d0*phi*alpha*(1.d0-beta*gamma/coshcplx(beta*x)**2)
     &                          /coshcplx(alpha*x)**2
     &        -2.d0*phi*gamma*beta**2*tanhcplx(alpha*x)*tanhcplx(beta*x)
     &                          /coshcplx(beta*x)**2
     &   +2.d0*alpha**2*phi*(x-gamma*tanhcplx(beta*x))*tanhcplx(alpha*x)
     &     /coshcplx(alpha*x)**2
     &  +(2.d0*x*beta**2*tanhcplx(beta*x)-2.d0*beta)/coshcplx(beta*x)**2

      RETURN
      END
!-----------------------------------------------------------------------
      FUNCTION dertrois(X, GAMMA, PHI,  ALPHA)
!
      IMPLICIT NONE
      COMPLEX(KIND(1.D0)) DERTROIS
      DOUBLE PRECISION, INTENT(IN) :: GAMMA, ALPHA
      COMPLEX(KIND(1.D0)), INTENT(IN) :: X, PHI
      COMPLEX(KIND(1.D0)) TANHCPLX,COSHCPLX
      EXTERNAL tanhcplx,coshcplx

      DOUBLE PRECISION BETA
      beta=1.d0-alpha

      dertrois=
     & -4.d0*phi*alpha*beta*gamma*(beta*tanhcplx(beta*x)
     &                              +alpha*tanhcplx(alpha*x))
     &                       /(coshcplx(alpha*x)**2*coshcplx(beta*x)**2)
     &       +4.d0*phi*alpha**2*tanhcplx(alpha*x)/(coshcplx(alpha*x)**2)
     &           -2.d0*phi*gamma*beta**2/(coshcplx(beta*x)**2)*
     &              ( alpha*tanhcplx(beta*x)/(coshcplx(alpha*x)**2)
     &                +beta*tanhcplx(alpha*x)/(coshcplx(beta*x)**2)
     &                -2.d0*beta*tanhcplx(alpha*x)*tanhcplx(beta*x)**2 )
     &           +2.d0*phi*alpha**2/(coshcplx(alpha*x)**2)*
     &       ( (1.d0-gamma*beta/(coshcplx(beta*x)**2))*tanhcplx(alpha*x)
     &                +(x-gamma*tanhcplx(beta*x))*
     &                 (alpha/(coshcplx(alpha*x)**2)
     &                  -2.d0*alpha*tanhcplx(alpha*x)**2) )
     &            +(2.d0*x*beta**3/(coshcplx(beta*x)**2)
     &               +6.d0*beta**2*tanhcplx(beta*x)
     &               -4.d0*x*beta**3*tanhcplx(beta*x)**2
     &              )/(coshcplx(beta*x)**2)
      RETURN
      END
!-----------------------------------------------------------------------
      FUNCTION tanhcplx(Z)
!
      USE declarations_special
      IMPLICIT NONE
      COMPLEX(KIND(1.D0)) TANHCPLX,III
      COMPLEX(KIND(1.D0)), INTENT(IN) :: Z

      iii = cmplx(0.d0,1.d0,kind=r8)

      tanhcplx= (exp(iii*z)+exp(-iii*z))/(exp(iii*z)-exp(-iii*z))
      RETURN
      END
!-----------------------------------------------------------------------
      FUNCTION coshcplx(Z)
!
      USE declarations_special
      IMPLICIT NONE
      COMPLEX(KIND(1.D0)) COSHCPLX,III
      COMPLEX(KIND(1.D0)), INTENT(IN) :: Z

      iii = cmplx(0.d0,1.d0,kind=r8)

      coshcplx= (exp(iii*z)+exp(-iii*z))*0.5d0
      RETURN
      END