breach.f

Go to the documentation of this file.

!                       ******************
                        SUBROUTINE breach
!                       ******************
!
!
!***********************************************************************
! TELEMAC2D   V6P2                                   03/08/2012
!***********************************************************************
!
!brief    MODIFICATION OF THE BOTTOM TOPOGRAPHY FOR BREACHES
!
!
!history  P. CHASSE (CETMEF) / C. COULET (ARTELIA)
!+        03/08/2012
!+        V6P2
!+        Creation
!
!history  Y.B. TADESSE (TUHH, INSTITUTE OF RIVER AND COASTAL ENGINEERING)
!+        14/02/2014
!+        V6P3R2
!+   Addition of later breach growth option
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_telemac2d
!
      USE declarations_special
      USE interface_parallel, ONLY : p_sum,p_max,p_min,
     &                               p_sum
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
      INTEGER I, J, K, N, M,END1, END2,CURNBR
      INTEGER ISTAT,VECZ
      INTEGER TEMPND(npoin)
      DOUBLE PRECISION ZC, ZW, ZB
      DOUBLE PRECISION AT1, AT2, AT3
!
      DOUBLE PRECISION X1, X2, Y1, Y2, DX, DY,DS1,DS2,END1X,END1Y
      DOUBLE PRECISION U1, U2, V1, V2, DELS,CURDIS,DIS,END2X,END2Y
      DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: XL, YL, XP, YP
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      IF (.NOT.dejalu_breach) THEN
        IF(debug.GT.0) WRITE(lu,*) 'CALLING LECBREACH'
        CALL lecbreach(t2d_files(t2dbrc)%LU)
        IF(debug.GT.0) WRITE(lu,*) 'BACK FROM LECBREACH'
        dejalu_breach=.true.
        WRITE (lu,*) 'READING BREACH DATA = OK'
!
        DO i = 1, nbrech
          zcrbr%R(i) = -huge(100.d0)
          DO j = 1, nbndbr%I(i)
            k = indbr%ADR(i)%P%I(j)
            IF(zf%R(k).GT.zcrbr%R(i)) THEN
              zcrbr%R(i) = zf%R(k)
            ENDIF
          ENDDO
          IF(ncsize.GT.1) THEN
            zcrbr%R(i) = p_max(zcrbr%R(i))
          ENDIF
        ENDDO
      ENDIF
!
      DO i = 1, nbrech
        IF (optero%I(i).EQ.1) THEN
          IF((optnbr%I(i).EQ.2).AND.(tdecbr%R(i).LT.0.d0)) THEN
            zc = 0.d0
            n = 0
            DO j = 1, nbndbr%I(i)
              k = indbr%ADR(i)%P%I(j)
              IF(h%R(k).GT.0.d0) THEN
                zw = zf%R(k)+h%R(k)
                zc = zc + zw
                n = n + 1
              ENDIF
            ENDDO
            IF(ncsize.GT.1) THEN
              n = p_sum(n)
              zc = p_sum(zc)
            ENDIF
            IF(n.GT.1) zc = zc/n
            IF(zc.GT.zdecbr%R(i)) THEN
              WRITE(lu,20) i, at
              tdecbr%R(i) = at
            ENDIF
          ENDIF
          IF((optnbr%I(i).EQ.3).AND.(tdecbr%R(i).LT.0.d0)) THEN
            IF(numpsd%I(i).GT.0) THEN
              zw = zf%R(numpsd%I(i)) + h%R(numpsd%I(i))
            ELSE
              zw = 0.d0
            ENDIF
!           CASE WHERE ONE OF THE ENDS IS NOT IN THE SUB-DOMAIN
            IF(ncsize.GT.1) THEN
              zw = p_max(zw)+p_min(zw)
            ENDIF
            IF(zw.GT.zdecbr%R(i)) THEN
              WRITE(lu,20) i, at
              tdecbr%R(i) = at
            ENDIF
          ENDIF
          at1 = tdecbr%R(i)
          at2 = at1 + durbr%R(i)
          IF(at1.GT.0.d0) THEN
            IF(at.GT.at1) THEN
              IF(at.GT.at2) THEN
                zb = zfinbr%R(i)
              ELSE
                zb = zcrbr%R(i)+(zfinbr%R(i)-zcrbr%R(i))/(at2-at1)
     &             *(at-at1)
              ENDIF
              DO j = 1, nbndbr%I(i)
                k = indbr%ADR(i)%P%I(j)
                zf%R(k)=min(zf%R(k), zb)
              ENDDO
            ENDIF
          ENDIF
        ELSEIF (optero%I(i).EQ.2) THEN
!       FIND NUMBER POINTS IN THE CURRENT BREACH REGION
          curdis=0.d0
          END1X=0.d0
          END1Y=0.d0
          END2X=0.d0
          END2Y=0.d0
          END1=0
          END2=0
          firstend: DO m=1,nponbr%I(i)-1
            curdis=curdis+pondsb%ADR(i)%P%R(m)
            IF (curdis .GT. ((finbrw%R(i)-curbrw%R(i))/2.d0)) THEN
              END1=M+1
              ds1= curdis - (finbrw%R(i)-curbrw%R(i))/2.d0
              dx = dkaxcr%ADR(i)%P%R(m+1) - dkaxcr%ADR(i)%P%R(m)
              dy = dkaycr%ADR(i)%P%R(m+1) - dkaycr%ADR(i)%P%R(m)
              u1 = dx/pondsb%ADR(i)%P%R(m)
              u2 = dy/pondsb%ADR(i)%P%R(m)
              IF(curbrw%R(i).LT.finbrw%R(i)) THEN
                END1X = DKAXCR%ADR(I)%P%R(M+1) - U1*ds1
                END1Y = DKAYCR%ADR(I)%P%R(M+1) - U2*ds1
              ELSEIF(curbrw%R(i).EQ.finbrw%R(i)) THEN
                END1X = DKAXCR%ADR(I)%P%R(1)- U1*
     &          (pondsb%ADR(i)%P%R(1)/10.d0)
                END1Y = DKAYCR%ADR(I)%P%R(1)- U2*
     &          (pondsb%ADR(i)%P%R(1)/10.d0)
              ELSE
                WRITE(lu,400)
                stop
              ENDIF
              dis=ds1
              secondend: DO j=m+1,nponbr%I(i)-1
                dis=dis+pondsb%ADR(i)%P%R(j)
                IF (dis .GE. curbrw%R(i)) THEN
                  END2=j
                  ds2=dis - curbrw%R(i)
                  dx = dkaxcr%ADR(i)%P%R(j+1) - dkaxcr%ADR(i)%P%R(j)
                  dy = dkaycr%ADR(i)%P%R(j+1) - dkaycr%ADR(i)%P%R(j)
                  u1 = dx/pondsb%ADR(i)%P%R(j)
                  u2 = dy/pondsb%ADR(i)%P%R(j)
                  IF(curbrw%R(i).LT.finbrw%R(i)) THEN
                    END2X = DKAXCR%ADR(I)%P%R(J+1) - U1*ds2
                    END2Y = DKAYCR%ADR(I)%P%R(J+1) - U2*ds2
                  ELSEIF (curbrw%R(i).EQ.finbrw%R(i)) THEN
                    END2X = DKAXCR%ADR(I)%P%R(NPONBR%I(I))+ U1*
     &              (pondsb%ADR(i)%P%R(nponbr%I(i)-1)/10.d0)
                    END2Y = DKAYCR%ADR(I)%P%R(NPONBR%I(I))+ U2*
     &              (pondsb%ADR(i)%P%R(nponbr%I(i)-1)/10.d0)
                  ELSE
                    WRITE(lu,400) i
                    stop
                  ENDIF
                  EXIT secondend
                ENDIF
              END DO secondend
!
              EXIT firstend
            ENDIF
          ENDDO firstend
400       FORMAT(1x, 'BREACH: BREACH WIDTH CANNOT BE GREATER ',
     &           'THAN FINAL WIDTH FOR BREACH: ', 1i6)
!         ALLOCATION OF LOCAL VARIABLE FOR BREACH DEFINITION
          istat = 0
          vecz=end2-end1+3
          ALLOCATE(xl(vecz), stat=istat)
          IF(istat.NE.0) THEN
            WRITE(lu,200) istat
            stop
          ENDIF
          ALLOCATE(yl(vecz), stat=istat)
          IF(istat.NE.0) THEN
            WRITE(lu,200) istat
            stop
          ENDIF
!
200       FORMAT(1x,'ERROR DURING ALLOCATION OF VECTOR: ',
     &           'ERROR CODE: ',1i6)
!
          xl(1)=end1x
          yl(1)=end1y
          xl(vecz)=end2x
          yl(vecz)=end2y
          DO j=2,vecz-1
            xl(j)=dkaxcr%ADR(i)%P%R(end1+j-2)
            yl(j)=dkaycr%ADR(i)%P%R(end1+j-2)
          ENDDO
!         SEARCH MESH POINTS INSIDE THE BREACH DOMAIN
          istat = 0
          ALLOCATE(xp(2*vecz), stat=istat)
          IF(istat.NE.0) THEN
            WRITE(lu,200) istat
            stop
          ENDIF
          ALLOCATE(yp(2*vecz), stat=istat)
          IF(istat.NE.0) THEN
            WRITE(lu,200) istat
            stop
          ENDIF
!
          x1 = xl(1)
          y1 = yl(1)
          x2 = xl(2)
          y2 = yl(2)
          dx = x2 - x1
          dy = y2 - y1
          dels=sqrt(dx*dx+dy*dy)
          IF(dels.GE.0.d0) THEN
            u1 = dx/dels
            u2 = dy/dels
          ELSE
            WRITE(lu,*) 'PROBLEM IN DEFINITION OF BREACH :',i
            CALL plante(1)
            stop
          ENDIF
          v1 = -u2
          v2 = u1
          xp(1)      = x1 + v1*polwdt%R(i)/2.d0
          yp(1)      = y1 + v2*polwdt%R(i)/2.d0
          xp(2*vecz) = x1 - v1*polwdt%R(i)/2.d0
          yp(2*vecz) = y1 - v2*polwdt%R(i)/2.d0
!
          DO m = 2,vecz
            x2 = xl(m)
            y2 = yl(m)
            dx = x2 - x1
            dy = y2 - y1
            dels=sqrt(dx*dx+dy*dy)
            IF(dels.GE.0.d0) THEN
              u1 = dx/dels
              u2 = dy/dels
            ELSE
              WRITE(lu,*) 'PROBLEM IN DEFINITION OF BREACH :',i
              CALL plante(1)
              stop
            ENDIF
            v1 = -u2
            v2 = u1
            xp(m)         = x2 + v1*polwdt%R(i)/2.d0
            yp(m)         = y2 + v2*polwdt%R(i)/2.d0
            xp(2*vecz-m+1) = x2 - v1*polwdt%R(i)/2.d0
            yp(2*vecz-m+1) = y2 - v2*polwdt%R(i)/2.d0
            x1=x2
            y1=y2
          ENDDO
          curnbr = 0
          DO m = 1, npoin
            IF(inpoly(mesh%X%R(m), mesh%Y%R(m), xp, yp, 2*vecz)) THEN
              curnbr = curnbr +1
              tempnd(curnbr) = m
            ENDIF
          ENDDO
!
          IF((optnbr%I(i).EQ.2).AND.(tdecbr%R(i).LT.0.d0)) THEN
            zc = 0.d0
            n = 0
            DO j = 1, curnbr
              k = tempnd(j)
              IF(h%R(k).GT.0.d0) THEN
                zw = zf%R(k)+h%R(k)
                zc = zc + zw
                n = n + 1
              ENDIF
            ENDDO
            IF(ncsize.GT.1) THEN
              n = p_sum(n)
              zc = p_sum(zc)
            ENDIF
            IF(n.GT.1) zc = zc/n
            IF(zc.GT.zdecbr%R(i)) THEN
              WRITE(lu,20) i, at
              tdecbr%R(i) = at
            ENDIF
          ENDIF
          IF((optnbr%I(i).EQ.3).AND.(tdecbr%R(i).LT.0.d0)) THEN
            IF(numpsd%I(i).GT.0) THEN
              zw = zf%R(numpsd%I(i)) + h%R(numpsd%I(i))
            ELSE
              zw = 0.d0
            ENDIF
!           CASE WHERE ONE OF THE ENDS IS NOT IN THE SUB-DOMAIN
            IF(ncsize.GT.1) THEN
              zw = p_max(zw)+p_min(zw)
            ENDIF
            IF(zw.GT.zdecbr%R(i)) THEN
              WRITE(lu,20) i, at
              tdecbr%R(i) = at
            ENDIF
          ENDIF
!
!         DURATION FOR THE VERTICAL EROSION IS TAKEN AS A TENTH
!         OF THE TOTAL BREACH DURATION
          at1 = tdecbr%R(i)
          at2 = at1 + durbr%R(i)
          at3= at1 + durbr%R(i)/10.d0
          IF(at1.GT.0.d0) THEN
            IF(at.GT.at1) THEN
              IF(at.GT.at3) THEN
                zb = zfinbr%R(i)
              ELSE
                zb = zcrbr%R(i)+(zfinbr%R(i)-zcrbr%R(i))/(at3-at1)
     &          *(at-at1)
              ENDIF
              IF(at.GT.at2) THEN
                curbrw%R(i)=finbrw%R(i)
              ELSE
                curbrw%R(i)=inibrw%R(i)+(finbrw%R(i)-inibrw%R(i))
     &           /(at2-at1)*(at-at1)
              ENDIF
              DO j = 1, curnbr
                k = tempnd(j)
                zf%R(k)=min(zf%R(k), zb)
              ENDDO
            ENDIF
          ENDIF
          DEALLOCATE(xl)
          DEALLOCATE(yl)
          DEALLOCATE(xp)
          DEALLOCATE(yp)
        ENDIF
      ENDDO
!
!-----------------------------------------------------------------------
!     MESSAGES
20    FORMAT(1x,'CREATION OF BREACH : ',i4,/,1x,
     &          'AT TIME : ',g16.7)
!-----------------------------------------------------------------------
!
      RETURN
      END