phbor.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE phbor
!                   ****************
!
!
!***********************************************************************
! ARTEMIS   V7P2                                     Nov 2016
!***********************************************************************
!
!brief    TRANSLATES THE BOUNDARY CONDITIONS SPECIFIED
!+                BY THE USER,
!+                I.E. COMPUTES THE COEFFICIENTS
!+                APHIR, APHII, ... FOR EACH BOUNDARY SEGMENT.
!
!history  J-M HERVOUET (LNH)
!+
!+
!+   LINKED TO BIEF 5.0
!
!history  D. AELBRECHT (LNH)
!+        21/08/2000
!+        V5P1
!+
!
!history  C. DENIS (SINETICS)
!+        18/03/2010
!+        V6P0
!+
!
!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  C.PEYRARD (LNHE)
!+        01/06/2012
!+        V6P0
!+   KSORT   : end of application to neighbours
!+   KINC    : end of application to the I+1 node
!+   SEGMENT : If a segment links a solid node to a liquid node,
!+             this segment is regarded as solid.
!+   BUG correction         : Full KINC boundaries taken into account
!+   NEW boundary condition : incident potential can be given by the user
!+   Parallel correction    : - End of HBT,CGT,CTT,KT,XT,YT tables
!+                              only   HB%R,CG%R, etc... are used.
!+                            - No use of NCSIZE variable.
!
!history  C.PEYRARD (LNHE)
!+        18/03/2014
!+        V7P0
!+   KSORT   : application to neighbours (best for automatic angles)
!
!history  N.DURAND (HRW)
!+        November 2016
!+        V7P2
!+   Phase calculation brought back inside PHBOR
!+   (automatic in the sense of applicable to general case, different
!+    from that using KPHREF)
!
!history  S.BOURBAN (HRW)
!+        23/01/2017
!+        V7P3
!+   Implementing a cumulative sorting of each boundary allowing for
!+   varying bathymetry along each boundary segment.
!
!history  N.DURAND (HRW)
!+        August 2017
!+        V7P3
!+   PI now defined in DECLARATIONS_ARTEMIS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_telemac
      USE declarations_artemis
!
      USE declarations_special
      IMPLICIT NONE
!
      INTEGER I,I0,IZ1,IZ2,IG,IGN,IGP,N
!
!
      DOUBLE PRECISION GRE, GIM
      DOUBLE PRECISION DDXGRE,DDYGRE,DDXGIM,DDYGIM
!

      DOUBLE PRECISION PHASEOI,X0,Y0
      DOUBLE PRECISION AUXI1,AUXIC,AUXIS
      DOUBLE PRECISION SDIST,MINDIST
!     DOUBLE PRECSION AUXI2
      INTRINSIC cos,sin
!
      DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-9
!
!-----------------------------------------------------------------------
!
      pi = 4.d0 * atan( 1.d0 )
      degrad = pi / 180.d0
      raddeg = 180.d0 / pi
!
!-----------------------------------------------------------------------
!
! INITIALISES LIDIR TO KSORT (A DIFFERENT VALUE FROM KENT)
! IN ORDER NOT TO TAKE NODES IMPOSED IN PRIDIH INTO ACCOUNT,
! WHEN IT HAS NOT BEEN REQUESTED.
!
      DO i=1,nptfr
        lidir%I(i) = ksort
!       BEWARE: IT IS ASSUMED HERE THAT NPTFRX=NPTFR
        lidir%I(i+nptfr) = ksort
        IF (lihbor%I(i).EQ.kent) THEN
          lihbor%I(i) = kinc
        ENDIF
        aphi1b%R(i) = 0.d0
        bphi1b%R(i) = 0.d0
        cphi1b%R(i) = 0.d0
        dphi1b%R(i) = 0.d0
        aphi2b%R(i) = 0.d0
        bphi2b%R(i) = 0.d0
        cphi2b%R(i) = 0.d0
        dphi2b%R(i) = 0.d0
        aphi3b%R(i) = 0.d0
        bphi3b%R(i) = 0.d0
        cphi3b%R(i) = 0.d0
        dphi3b%R(i) = 0.d0
        aphi4b%R(i) = 0.d0
        bphi4b%R(i) = 0.d0
        cphi4b%R(i) = 0.d0
        dphi4b%R(i) = 0.d0
        cgrx1b%R(i) = 0.d0
        cgry1b%R(i) = 0.d0
        dgrx1b%R(i) = 0.d0
        dgry1b%R(i) = 0.d0
      ENDDO
!
!       -------------------------------------------------
!             COEFFICIENTS FOR A INCIDENT WAVE
!       -------------------------------------------------
!
!-----------------------------------------------------------------------
!
!     PROCESSING EACH BOUNDARY SEGMENT INDEPEDENTLY,
!     STARTING WITH THE FIRST NODE HIT BY THE INCOMING WAVE
!
      DO n = 1,nfrliq
!
!     ***********************************************************
!     LOOKING FOR AN APPROPRIATE STARTING POINT: I0 (AND END POINTS IZS)
!
        i0 = 0
        DO i=1,nptfr
          IF( numliq%I(i).EQ.n ) THEN
!     SIGNED DISTANCE BETWEEN THE CREST OF THE WAVE AND THE BOUNDARY
!     ( TETAB%R(I) DEFINED AS WAVE PROPAGATING TOWARDS )
            ig   = mesh%NBOR%I(i)
            sdist = cos(tetab%R(i)*degrad)*( x_phref*y_phref-x(ig) )
     &            + sin(tetab%R(i)*degrad)*( x_phref*y_phref-y(ig) )
            IF( i0.EQ.0 ) THEN
              i0 = i
              mindist = sdist
            ELSE
              IF( sdist.GT.mindist ) THEN
                i0 = i
                mindist = sdist
              ENDIF
            ENDIF
          ENDIF
        ENDDO
!
        IF( i0.EQ.0 ) cycle
!
        iz1 = mesh%KP1BOR%I(i0)
        iz2 = mesh%KP1BOR%I(i0+nptfr)
        DO WHILE( numliq%I(iz1).EQ.n .AND. numliq%I(iz2).EQ.n )
          iz1 = mesh%KP1BOR%I(iz1)
          iz2 = mesh%KP1BOR%I(iz2+nptfr)
          IF( iz1.EQ.iz2 ) EXIT
        ENDDO
        IF( numliq%I(iz1).NE.n ) THEN
          DO WHILE( numliq%I(iz2).EQ.n )
            iz2 = mesh%KP1BOR%I(iz2+nptfr)
          ENDDO
        ENDIF
        IF( numliq%I(iz2).NE.n ) THEN
          DO WHILE( numliq%I(iz1).EQ.n )
            iz1 = mesh%KP1BOR%I(iz1)
          ENDDO
        ENDIF
!
!       CASE OF A BOUNDARY SPLIT OVER MORE THAN ONE PROCESSOR
! /!\   THIS SHOULD NOT HAPPEN AGAIN BECAUSE ARTEMIS'S SPECIAL TREATMENT
!        IF( NCSIZE.GT.1 ) THEN
!          SDIST = MDIST
!          MDIST = P_MIN( MDIST )
!          IG = 0
!          IF( ABS(SDIST-MDIST).LT.CHOUIA ) IG = I0
!          I0 = P_MAX( I0 )
!        ENDIF
!
!     ***********************************************************
!     DUAL WAVE PROPAGATION ALONG THE BOUNDARY, STARTING FROM I0
!
        IF (lihbor%I(i0).EQ.kinc) THEN
!
          WRITE(lu,112) n,mesh%NBOR%I(i0),
     &      x(mesh%NBOR%I(i0)),y(mesh%NBOR%I(i0))
 112      FORMAT(/,1x,'PHBOR: BOUNDARY ',i3,
     &      ' IS INCIDENT - FIRST IMPACT NODE ',i6,
     &      '[',f12.4,',',f12.4,']')
!
!     1.  FOLLOWING THE BOUNDARY NODES GOING ONE WAY (INCL. I0)
          i  = i0
          ig = mesh%NBOR%I(i)
          x0 = x_phref
          y0 = y_phref
          phaseoi = 0.d0
          DO WHILE( numliq%I(i).EQ.n )
!           -- LOCAL CONSTANTS
            auxic = cos(tetab%R(i)*degrad)
            auxis = sin(tetab%R(i)*degrad)
            auxi1 = grav/omega * hb%R(i)/2.d0
!           -- AUTOMATIC PHASE CALCULATION (TETAB HAVE TO BE UNIFORM ON THE BOUNDARY)
            IF(lphaseauto) THEN
              phaseoi =
     &          kphref*auxic*( x(ig)-x0 ) + kphref*auxis*( y(ig)-y0 )
            ELSE
              phaseoi = phaseoi +
     &          k%R(ig)*auxic*( x(ig)-x0 ) + k%R(ig)*auxis*( y(ig)-y0 )
!             -- INCREMENT WITHIN
              x0=x(ig)
              y0=y(ig)
            ENDIF
!           -- INCIDENT WAVE --> INCIDENT POTENTIAL (REAL, IMAGINAR)
            gre= auxi1*sin( alfap%R(i)*degrad + phaseoi)
            gim=-auxi1*cos( alfap%R(i)*degrad + phaseoi)
!           -- INCIDENT WAVE --> GRADIENTS          (REAL, IMAGINAR)
            ddxgre= auxi1*cos(alfap%R(i)*degrad+phaseoi)*auxic*k%R(ig)
            ddygre= auxi1*cos(alfap%R(i)*degrad+phaseoi)*auxis*k%R(ig)
            ddxgim= auxi1*sin(alfap%R(i)*degrad+phaseoi)*auxic*k%R(ig)
            ddygim= auxi1*sin(alfap%R(i)*degrad+phaseoi)*auxis*k%R(ig)
!           -- MATRIX AM AND BM COEFFICIENTS
            aphi1b%R(i) = - k%R(ig) * c%R(ig) * cg%R(ig)
     &                   * cos(tetap%R(i)*degrad)
            bphi1b%R(i) = 0.d0
!           -- SECOND MEMBER CV1 AND CV2 COEFFICIENTS
!           i * K * Gamma (multiplied by "- cos THETAP" in BERKHO)
            cphi1b%R(i)  = ( -gim*k%R(ig) ) *c%R(ig)*cg%R(ig)
            dphi1b%R(i)  = (  gre*k%R(ig) ) *c%R(ig)*cg%R(ig)
!           -- GRAD(Gamma) will be used in BERKHO ...
            cgrx1b%R(i)=   ( ddxgre ) *c%R(ig)*cg%R(ig)
            cgry1b%R(i)=   ( ddygre ) *c%R(ig)*cg%R(ig)
            dgrx1b%R(i)=   ( ddxgim ) *c%R(ig)*cg%R(ig)
            dgry1b%R(i)=   ( ddygim ) *c%R(ig)*cg%R(ig)
!           -- WHILE INCREMENT
            IF( i.EQ.iz1 ) EXIT
            i  = mesh%KP1BOR%I(i)
            ig = mesh%NBOR%I(i)
          ENDDO
!
!     2.  FOLLOWING THE BOUNDARY GOING THE OTHER WAY (INCL. I0)
          i  = i0
          ig = mesh%NBOR%I(i)
          x0 = x_phref
          y0 = y_phref
          phaseoi = 0.d0
          DO WHILE( numliq%I(i).EQ.n )
!           -- LOCAL CONSTANTS
            auxic = cos(tetab%R(i)*degrad)
            auxis = sin(tetab%R(i)*degrad)
            auxi1 = grav/omega * hb%R(i)/2.d0
!           -- AUTOMATIC PHASE CALCULATION (TETAB HAVE TO BE UNIFORM ON THE BOUNDARY)
            IF(lphaseauto) THEN
              phaseoi =
     &          kphref*auxic*( x(ig)-x0 ) + kphref*auxis*( y(ig)-y0 )
            ELSE
              phaseoi = phaseoi +
     &          k%R(ig)*auxic*( x(ig)-x0 ) + k%R(ig)*auxis*( y(ig)-y0 )
!             -- INCREMENT WITHIN
              x0=x(ig)
              y0=y(ig)
            ENDIF
!           -- INCIDENT WAVE --> INCIDENT POTENTIAL (REAL, IMAGINAR)
            gre= auxi1*sin( alfap%R(i)*degrad + phaseoi)
            gim=-auxi1*cos( alfap%R(i)*degrad + phaseoi)
!           -- INCIDENT WAVE --> GRADIENTS          (REAL, IMAGINAR)
            ddxgre= auxi1*cos(alfap%R(i)*degrad+phaseoi)*auxic*k%R(ig)
            ddygre= auxi1*cos(alfap%R(i)*degrad+phaseoi)*auxis*k%R(ig)
            ddxgim= auxi1*sin(alfap%R(i)*degrad+phaseoi)*auxic*k%R(ig)
            ddygim= auxi1*sin(alfap%R(i)*degrad+phaseoi)*auxis*k%R(ig)
!           -- MATRIX AM AND BM COEFFICIENTS
            aphi1b%R(i) = - k%R(ig) * c%R(ig) * cg%R(ig)
     &                   * cos(tetap%R(i)*degrad)
            bphi1b%R(i) = 0.d0
!           -- SECOND MEMBER CV1 AND CV2 COEFFICIENTS
!           i * K * Gamma (multiplied by "- cos THETAP" in BERKHO)
            cphi1b%R(i)  = ( -gim*k%R(ig) ) *c%R(ig)*cg%R(ig)
            dphi1b%R(i)  = (  gre*k%R(ig) ) *c%R(ig)*cg%R(ig)
!           -- GRAD(Gamma) will be used in BERKHO ...
            cgrx1b%R(i)=   ( ddxgre ) *c%R(ig)*cg%R(ig)
            cgry1b%R(i)=   ( ddygre ) *c%R(ig)*cg%R(ig)
            dgrx1b%R(i)=   ( ddxgim ) *c%R(ig)*cg%R(ig)
            dgry1b%R(i)=   ( ddygim ) *c%R(ig)*cg%R(ig)
!           -- WHILE INCREMENT
            IF( i.EQ.iz2 ) EXIT
            i  = mesh%KP1BOR%I(i+nptfr)
            ig = mesh%NBOR%I(i)
          ENDDO

        ENDIF
!
!     ***********************************************************
!
      ENDDO
!
!     INITIALISATION OF PHASE VARIABLES FOR AUTOMATIC CALCULATION
      x0      = x_phref
      y0      = y_phref
      phaseoi = 0d0
!
      DO i=1,nptfr
!
!       ********************************
!       GLOBAL NUMBER OF THE BOUNDARY NODE I
!       ********************************
!
        ig  = mesh%NBOR%I(i)
!
!       ******************************************
!       GLOBAL NUMBER OF THE BOUNDARY NODE PRECEDING I
!       ******************************************
!
        ign = mesh%NBOR%I(mesh%KP1BOR%I(i+nptfr))
!
!       ****************************************
!       GLOBAL NUMBER OF THE BOUNDARY NODE FOLLOWING I
!       ****************************************
!
        igp = mesh%NBOR%I(mesh%KP1BOR%I(i))

!       -------------------------------------------------
!             COEFFICIENTS FOR AN INCIDENT WAVE
!       -------------------------------------------------
!             ALREADY TREATED ABOVE
!
!        IF (LIHBOR%I(I).EQ.KINC) THEN
!
!        ENDIF
!
!        -------------------------------------------------
!             COEFFICIENTS FOR AN INCIDENT POTENTIAL
!        -------------------------------------------------
!
        IF (lihbor%I(i).EQ.kpot) THEN
!
!------------ POTENTIAL (REAL, IMAGINAR)
          gre=prb%R(i)
          gim=pib%R(i)
! ----------- GRADIENTS (REAL, IMAGINAR)
          ddxgre= ddxprb%R(i)
          ddygre= ddyprb%R(i)
          ddxgim= ddxpib%R(i)
          ddygim= ddypib%R(i)
! --------COEFFICIENTS
! --      MATRIX AM AND BM COEFFICIENTS
          aphi1b%R(i) = - k%R(ig) * c%R(ig) * cg%R(ig)
     &                  * cos(tetap%R(i)*degrad)
          bphi1b%R(i) = 0.d0
! --      SECOND MEMBER CV1 AND CV2 COEFFICIENTS
! ----    i * K * Gamma (multiplied by "- cos THETAP" in BERKHO)
          cphi1b%R(i)  = ( -gim*k%R(ig) ) *c%R(ig)*cg%R(ig)
          dphi1b%R(i)  = (  gre*k%R(ig) ) *c%R(ig)*cg%R(ig)
! ----    GRAD(Gamma) will be used in BERKHO...
          cgrx1b%R(i)=   ( ddxgre ) *c%R(ig)*cg%R(ig)
          cgry1b%R(i)=   ( ddygre ) *c%R(ig)*cg%R(ig)
          dgrx1b%R(i)=   ( ddxgim ) *c%R(ig)*cg%R(ig)
          dgry1b%R(i)=   ( ddygim ) *c%R(ig)*cg%R(ig)
        ENDIF
!         -------------------------------------------------
!         COEFFICIENTS FOR A FREE EXIT BOUNDARY SEGMENT
!         -------------------------------------------------
!
        IF(lihbor%I(i).EQ.ksort) THEN
!
          aphi2b%R(i)  = - k%R(ig) * c%R(ig) * cg%R(ig)
     &                   * cos(tetap%R(i)*degrad)
          bphi2b%R(i)  = 0.d0
          cphi2b%R(i)  = 0.d0
          dphi2b%R(i)  = 0.d0
!
        ELSEIF (lihbor%I(mesh%KP1BOR%I(i)).EQ.ksort) THEN
          aphi2b%R(i)  = - k%R(ig) * c%R(ig) * cg%R(ig)
     &                   * cos(tetap%R(mesh%KP1BOR%I(i))*degrad)
          bphi2b%R(i)  = 0.d0
          cphi2b%R(i)  = 0.d0
          dphi2b%R(i)  = 0.d0
!
        ELSEIF (lihbor%I(mesh%KP1BOR%I(i+nptfr)).EQ.ksort) THEN
          aphi2b%R(i)  = - k%R(ig) * c%R(ig) * cg%R(ig)
     &                  * cos(tetap%R(mesh%KP1BOR%I(i+nptfr))*degrad)
          bphi2b%R(i)  = 0.d0
          cphi2b%R(i)  = 0.d0
          dphi2b%R(i)  = 0.d0
!
        ELSE
          aphi2b%R(i)  = 0.d0
          bphi2b%R(i)  = 0.d0
          cphi2b%R(i)  = 0.d0
          dphi2b%R(i)  = 0.d0
        ENDIF
!
!       -------------------------------------------
!       COEFFICIENTS FOR A SOLID BOUNDARY SEGMENT
!       -------------------------------------------
!
        IF (lihbor%I(i).EQ.klog) THEN
!
          auxi1 = k%R(ig) * c%R(ig) * cg%R(ig) *
     &             cos(tetap%R(i)*degrad) /
     &             (1.d0 + rp%R(i)*rp%R(i) +
     &             2.d0*rp%R(i)*cos(alfap%R(i)*degrad))
!
          aphi3b%R(i) = - (1.d0 - rp%R(i) * rp%R(i) ) * auxi1
!
          bphi3b%R(i) = 2.d0*rp%R(i)*sin(alfap%R(i)*degrad) * auxi1
!
          cphi3b%R(i)  = 0.d0
!
          dphi3b%R(i)  = 0.d0
!
        ELSEIF (lihbor%I(mesh%KP1BOR%I(i)).EQ.klog) THEN
          auxi1 = k%R(ig) * c%R(ig) * cg%R(ig) *
     &      cos(tetap%R(mesh%KP1BOR%I(i))*degrad) /
     &      (1.d0 + rp%R(mesh%KP1BOR%I(i))*rp%R(mesh%KP1BOR%I(i))
     &      +2.d0 * rp%R(mesh%KP1BOR%I(i))*
     &       cos(alfap%R(mesh%KP1BOR%I(i))*degrad))
!
          aphi3b%R(i) = - (1.d0-rp%R(mesh%KP1BOR%I(i))*
     &      rp%R(mesh%KP1BOR%I(i))) * auxi1
!
          bphi3b%R(i) = 2.d0*rp%R(mesh%KP1BOR%I(i))
     &                * sin(alfap%R(mesh%KP1BOR%I(i))*degrad) * auxi1
!
          cphi3b%R(i)  = 0.d0
!
          dphi3b%R(i)  = 0.d0
!
        ELSEIF (lihbor%I(mesh%KP1BOR%I(i+nptfr)).EQ.klog) THEN
          auxi1 = k%R(ig) * c%R(ig) * cg%R(ig) *
     &     cos(tetap%R(mesh%KP1BOR%I(i+nptfr))*degrad) /
     &     (1.d0 + rp%R(mesh%KP1BOR%I(i+nptfr))
     &      *rp%R(mesh%KP1BOR%I(i+nptfr))
     &      +2.d0 * rp%R(mesh%KP1BOR%I(i+nptfr))*
     &       cos(alfap%R(mesh%KP1BOR%I(i+nptfr))*degrad))
!
          aphi3b%R(i) = - (1.d0-rp%R(mesh%KP1BOR%I(i+nptfr))*
     &       rp%R(mesh%KP1BOR%I(i+nptfr))) * auxi1
!
          bphi3b%R(i) = 2.d0*rp%R(mesh%KP1BOR%I(i+nptfr))
     &      * sin(alfap%R(mesh%KP1BOR%I(i+nptfr))*degrad) * auxi1
!
          cphi3b%R(i)  = 0.d0
!
          dphi3b%R(i)  = 0.d0
!
        ELSE
          aphi3b%R(i)  = 0.d0
!
          bphi3b%R(i)  = 0.d0
!
          cphi3b%R(i)  = 0.d0
!
          dphi3b%R(i)  = 0.d0
!
        ENDIF
!
!        -------------------------------------------------
!        COEFFICIENTS FOR AN IMPOSED WAVE BOUNDARY SEGMENT
!        -------------------------------------------------
!DA      -----------------------------------
!DA      KEPT FOR MEMORY !
!DA      -----------------------------------
!DA
!DA         IF (LIHBOR(I).EQ.KENT) THEN
!DA         AUXIC      = COS(TETAB(I)*DEGRAD)
!DA         AUXIS      = SIN(TETAB(I)*DEGRAD)
!DA         AUXI1      = GRAV/OMEGA * HB(I)/2.D0 * C(IG) * CG(IG) *
!DA     *                K(IG) * ( AUXIC *XSGBOR(I) +
!DA     *                          AUXIS *YSGBOR(I) )
!DA         AUXI2      = K(IG) * ( X(IG)*AUXIC +
!DA     *                          Y(IG)*AUXIS )
!DA
!DA         APHI4B(I)  = 0.D0
!DA
!DA         BPHI4B(I)  = 0.D0
!DA
!DA         CPHI4B(I)  = AUXI1 * COS( AUXI2 )
!DA
!DA         DPHI4B(I)  = AUXI1 * SIN( AUXI2 )
!DA
!DA       VALUES IMPOSED AT THE NODES OF A KENT SEGMENT
!DA         LIDIR(I)         = KENT
!DA
!DA         AUXI1 = GRAV/OMEGA * HB(I)/2.D0
!DA         AUXI2 = K(IG) * (X(IG)*AUXIC +
!DA     *                    Y(IG)*AUXIS )
!DA
!DA            PHIRB(I) =   AUXI1 * SIN( AUXI2 )
!DA            PHIIB(I) = - AUXI1 * COS( AUXI2 )
!DA         ENDIF
!
!
      ENDDO
!-----------------------------------------------------------------------
!
!
      RETURN
      END SUBROUTINE