maxslope.f

Go to the documentation of this file.

!                   *******************
                    SUBROUTINE maxslope
!                   *******************
!
     &(slope,zf,zr,xel,yel,nelem,nelmax,npoin,ikle,evol,unsv2d,mesh,
     & zfcl_ms,avail,nomblay,nsicla)
!
!***********************************************************************
! SISYPHE   V7P2                                   14/11/2016
!***********************************************************************
!
!brief    COLLAPSE OF SAND WITH A SLOPE GREATER THAN A
!+                STABILITY CRITERION.
!+        For more explanation see release notes 5.8
!
!history  J-M HERVOUET (LNH)
!+        16/11/2007
!+        V5P8
!+
!
!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  J-M HERVOUET (EDF R&D, LNHE)
!+        08/03/2013
!+        V6P3
!+   Now possible with several classes of sediment.
!
!history  MICHIEL KNAAPEN (HRW)
!+        08/12/2016
!+        V7P2
!+        Improvements of the implementation of the morphological factor
!+        in the subroutine. Solution of stability issues. For futher
!+        details see the documentation.
!+
!history R KOPMANN (BAW)
!+        18/05/2018
!+        V7P3
!+   Angle will reduced to avoid changes below active layer only in case
!+   of multi grain
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| EVOL           |<->| WORK ARRAY, THEN EVOLUTION DUE TO SLIDE
!| IKLE           |-->| CONNECTIVITY TABLE
!| MESH           |-->| MESH STRUCTURE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS IN THE MESH
!| SLOPE          |-->| MAXIMUM SLOPE IN DEGREES
!| UNSV2D         |-->| INVERSE OF INTEGRAL OF BASES
!| XEL,YEL        |-->| MESH COORDINATES PER ELEMENT
!| ZF             |<->| BOTTOM THAT WILL BE MODIFIED
!| ZR             |-->| NON ERODABLE BED
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_sisyphe, ONLY : mofac
!
      USE interface_sisyphe, ex_maxslope => maxslope
!
      USE declarations_sisyphe, ONLY : es
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN,NOMBLAY,NSICLA
      INTEGER, INTENT(IN) :: IKLE(nelmax,3)
!
      DOUBLE PRECISION, INTENT(IN   ) :: SLOPE
      DOUBLE PRECISION, INTENT(INOUT) :: ZF(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: ZR(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: XEL(nelmax,3),YEL(nelmax,3)
      DOUBLE PRECISION, INTENT(INOUT) :: AVAIL(npoin,nomblay,nsicla)
!
      TYPE(bief_obj), INTENT(INOUT)   :: EVOL,ZFCL_MS
      TYPE(bief_obj), INTENT(IN)      :: UNSV2D
      TYPE(bief_mesh) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,I1,I2,I3,I,IG1,IG2,IR1,IR2,J
      DOUBLE PRECISION X2,X3,Y2,Y3,Z2,Z3,A,B,L,ZC,DEUXSURF,TANSL
      DOUBLE PRECISION Q(3),QG1,QG2,QR1,QR2
      DOUBLE PRECISION PI
      DOUBLE PRECISION EZ1,EZ2,EZ3
!
      LOGICAL CASE2
!
      INTRINSIC sqrt,min,max,tan
!
!-----------------------------------------------------------------------
!
      pi = 4.d0 * atan( 1.d0 )
      tansl = tan( pi*slope/180.d0 )
!
!     INITIALISES THE RIGHT-HAND SIDE EVOL TO ZERO
!
      CALL cpstvc(unsv2d,evol)
      CALL os('X=0     ',x=evol)
!
!     ONE CLASS VERSION
!
      IF(nsicla.EQ.1) THEN
!
!       LOOP ON ELEMENTS
!
        DO ielem=1,nelem
!
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
!
          x2=xel(ielem,2)
          x3=xel(ielem,3)
          y2=yel(ielem,2)
          y3=yel(ielem,3)
          z2=zf(i2)-zf(i1)
          z3=zf(i3)-zf(i1)
!
!         TWICE THE TRIANGLE AREA
!
          deuxsurf=x2*y3-x3*y2
!
!         AVERAGE BOTTOM IN THE ELEMENT
!
          zc=(zf(i1)+zf(i2)+zf(i3))/3.d0
!
!         COMPONENTS OF BOTTOM GRADIENT
!
          a=(z2*y3-z3*y2)/deuxsurf
          b=(z3*x2-z2*x3)/deuxsurf
!
!         CORRECTING FACTOR ON SLOPE
!
          l=min(1.d0,tansl/max(sqrt(a**2+b**2),1.d-8))
!
!         L LIMITED DUE TO NON-ERODABLE BEDS : ZF MUST NOT GO BELOW ZR
!
          IF(zf(i1).GT.zc) l=max(l,(zr(i1)-zc)/max(zf(i1)-zc,1.d-8))
          IF(zf(i2).GT.zc) l=max(l,(zr(i2)-zc)/max(zf(i2)-zc,1.d-8))
          IF(zf(i3).GT.zc) l=max(l,(zr(i3)-zc)/max(zf(i3)-zc,1.d-8))
!
!         BUILDS THE RIGHT-HAND SIDE
!
!         HERE THE EVOLUTIONS ARE MULTIPLIED BY SURFAC/3
!         BECAUSE THE REAL EVOLUTION TAKING INTO ACCOUNT OTHER ELEMENTS
!         WILL NEED A FACTOR (SURFAC/3)/(INTEGRAL OF BASIS)
!
          evol%R(i1)=evol%R(i1)+(1.d0-l)*(zc-zf(i1))*deuxsurf/6.d0
          evol%R(i2)=evol%R(i2)+(1.d0-l)*(zc-zf(i2))*deuxsurf/6.d0
          evol%R(i3)=evol%R(i3)+(1.d0-l)*(zc-zf(i3))*deuxsurf/6.d0
!
        ENDDO
!
      ELSE
!
!       MULTI-CLASS VERSION
!
!       INITIALING TO 0. THE EVOLUTIONS DUE TO SLIDE FOR EACH CLASS
!
        DO i=1,nsicla
          CALL os('X=0     ',x=zfcl_ms%ADR(i)%P)
        ENDDO
!
!       LOOP ON ELEMENTS
!
        DO ielem=1,nelem
!
          i1=ikle(ielem,1)
          i2=ikle(ielem,2)
          i3=ikle(ielem,3)
!
          x2=xel(ielem,2)
          x3=xel(ielem,3)
          y2=yel(ielem,2)
          y3=yel(ielem,3)
          z2=zf(i2)-zf(i1)
          z3=zf(i3)-zf(i1)
!
!         TWICE THE TRIANGLE AREA
!
          deuxsurf=x2*y3-x3*y2
!
!         AVERAGE BOTTOM IN THE ELEMENT
!
          zc=(zf(i1)+zf(i2)+zf(i3))/3.d0
!
!         COMPONENTS OF BOTTOM GRADIENT
!
          a=(z2*y3-z3*y2)/deuxsurf
          b=(z3*x2-z2*x3)/deuxsurf
!
!         CORRECTING FACTOR ON SLOPE
!
          l=min(1.d0,tansl/max(sqrt(a**2+b**2),1.d-8))
!
!         L LIMITED DUE TO NON-ERODABLE BEDS AND ACTIVE LAYER THICKNESS:
!         ZF MUST NOT GO BELOW ZR OR ACTIVE LAYER
!
!          IF(ZF(I1).GT.ZC) L=MAX(L,(ZR(I1)-ZC)/MAX(ZF(I1)-ZC,1.D-8))
!          IF(ZF(I2).GT.ZC) L=MAX(L,(ZR(I2)-ZC)/MAX(ZF(I2)-ZC,1.D-8))
!          IF(ZF(I3).GT.ZC) L=MAX(L,(ZR(I3)-ZC)/MAX(ZF(I3)-ZC,1.D-8))
          ez1 = max(zr(i1),zf(i1)-es(i1,1))
          ez2 = max(zr(i2),zf(i2)-es(i2,1))
          ez3 = max(zr(i3),zf(i3)-es(i3,1))
          IF(zf(i1).GT.zc) l=max(l,(ez1-zc)/max(zf(i1)-zc,1.d-8))
          IF(zf(i2).GT.zc) l=max(l,(ez2-zc)/max(zf(i2)-zc,1.d-8))
          IF(zf(i3).GT.zc) l=max(l,(ez3-zc)/max(zf(i3)-zc,1.d-8))
!
!         BUILDS THE RIGHT-HAND SIDE
!
!         HERE THE EVOLUTIONS ARE MULTIPLIED BY SURFAC/3
!         BECAUSE THE REAL EVOLUTION TAKING INTO ACCOUNT OTHER ELEMENTS
!         WILL NEED A FACTOR (SURFAC/3)/(INTEGRAL OF BASIS)
!
!         FIRST IN TERMS OF QUANTITIES BROUGHT TO POINTS
!
          q(1)=(1.d0-l)*(zc-zf(i1))*deuxsurf/6.d0
          q(2)=(1.d0-l)*(zc-zf(i2))*deuxsurf/6.d0
          q(3)=(1.d0-l)*(zc-zf(i3))*deuxsurf/6.d0
!
          evol%R(i1)=evol%R(i1)+q(1)
          evol%R(i2)=evol%R(i2)+q(2)
          evol%R(i3)=evol%R(i3)+q(3)
!
!         TAKING INTO ACCOUNT THE QUANTITIES TO UPDATE ZFCL_MS
!         IG1 AND IG2 : POINTS THAT GIVE
!         IR1 AND IR2 : POINTS THAT RECEIVE
!         CASE2: TWO POINTS GIVE TO THE THIRD ONE (THE OTHER CASE IS
!                ONE POINT GIVES TO THE TWO OTHERS)
          case2=.false.
!
!         PARAMETERISING TO REDUCE THE 6 CASES TO 2
!
          IF(q(1).GE.0.d0) THEN
            IF(q(2).GE.0.d0) THEN
!             3 GIVES TO 1 AND 2
              ig1=i3
              qg1=q(3)
              ir1=i1
              qr1=q(1)
              ir2=i2
              qr2=q(2)
            ELSE
              IF(q(3).GE.0.d0) THEN
!               2 GIVES TO 1 AND 3
                ig1=i2
                qg1=q(2)
                ir1=i1
                qr1=q(1)
                ir2=i3
                qr2=q(3)
              ELSE
!               2 AND 3 GIVE TO 1
                ig1=i2
                qg1=q(2)
                ig2=i3
                qg2=q(3)
                ir1=i1
                qr1=q(1)
                case2=.true.
              ENDIF
            ENDIF
          ELSE
            IF(q(2).GT.0.d0) THEN
              IF(q(3).GT.0.d0) THEN
!               1 GIVES TO 2 AND 3
                ig1=i1
                qg1=q(1)
                ir1=i2
                qr1=q(2)
                ir2=i3
                qr2=q(3)
              ELSE
!               1 AND 3 GIVE TO 2
                ig1=i1
                qg1=q(1)
                ig2=i3
                qg2=q(3)
                ir1=i2
                qr1=q(2)
                case2=.true.
              ENDIF
            ELSE
!             1 AND 2 GIVE TO 3
              ig1=i1
              qg1=q(1)
              ig2=i2
              qg2=q(2)
              ir1=i3
              qr1=q(3)
              case2=.true.
            ENDIF
          ENDIF
!
          IF(case2) THEN
!
!           THE TWO DONNORS CASE : IG1 AND IG2 GIVE TO IR1
!           ZFCL_MS IS HERE VOLUMES
!
            DO i=1,nsicla
              zfcl_ms%ADR(i)%P%R(ig1)=zfcl_ms%ADR(i)%P%R(ig1)
     &                               +qg1*avail(ig1,1,i)
              zfcl_ms%ADR(i)%P%R(ig2)=zfcl_ms%ADR(i)%P%R(ig2)
     &                               +qg2*avail(ig2,1,i)
              zfcl_ms%ADR(i)%P%R(ir1)=zfcl_ms%ADR(i)%P%R(ir1)
     &                               -qg1*avail(ig1,1,i)
     &                               -qg2*avail(ig2,1,i)
            ENDDO
!
          ELSE
!
!           THE ONE DONNOR CASE : IG1 GIVES TO IR1 AND IR2
!           ZFCL_MS IS HERE VOLUMES
!
            DO i=1,nsicla
              zfcl_ms%ADR(i)%P%R(ig1)=zfcl_ms%ADR(i)%P%R(ig1)
     &                               +qg1*avail(ig1,1,i)
              zfcl_ms%ADR(i)%P%R(ir1)=zfcl_ms%ADR(i)%P%R(ir1)
     &                               +qr1*avail(ig1,1,i)
              zfcl_ms%ADR(i)%P%R(ir2)=zfcl_ms%ADR(i)%P%R(ir2)
     &                               +qr2*avail(ig1,1,i)
            ENDDO
!
          ENDIF
!
        ENDDO
!
!       ADDING VOLUMES IN PARALLEL
!
        IF(ncsize.GT.1) THEN
          DO i=1,nsicla
            CALL parcom(zfcl_ms%ADR(i)%P,2,mesh)
          ENDDO
        ENDIF
!
!       FINAL DIVISION BY THE INTEGRAL OF BASES: VOLUMES CHANGED INTO
!       BED VARIATIONS (LIKE EVOL BELOW)
!
        DO i=1,nsicla
          DO j=1,npoin
            zfcl_ms%ADR(i)%P%R(j)=zfcl_ms%ADR(i)%P%R(j)*
     &                unsv2d%R(j)/mofac/10.d0
          ENDDO
        ENDDO
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!     FINAL RESOLUTION
!
      IF(ncsize.GT.1) THEN
        CALL parcom(evol,2,mesh)
      ENDIF
!
!     FINAL DIVISION BY THE INTEGRAL OF BASES: QUANTITIES CHANGED INTO
!     ELEVATIONS
!
      DO i=1,npoin
        evol%R(i)=evol%R(i)*unsv2d%R(i)/mofac/10.d0
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END