layer.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE layer
!                   ****************
!
     &(zfcl_w,nlayer,zr,zf,estrat,elay,masbas,acladm,nsicla,npoin,
     & elay0,voltot,es,avail,const_alayer,estratnew,nlaynew)
!
!***********************************************************************
! SISYPHE   V6P2                                   21/07/2011
!***********************************************************************
!
!brief    COMPUTES AVAIL FOR EACH CLASS AND EACH LAYER;
!+                NEW STRATUM THICKNESS ESTRAT.
!+
!+            ACTIVE LAYER IS LAYER 1, IT IS KEPT AT A PRESCRIBED
!+                HEIGHT OF ELAY0  PROVIDED IT IS POSSIBLE
!+
!+            STRATUM IS LAYER 2 OF HEIGHT ESTRAT.
!
!history  MATTHIEU GONZALES DE LINARES
!+        2002
!+
!+
!
!history  JMH
!+        16/09/2009
!+
!+   AVAIL(NPOIN,10,NSICLA)
!
!history  JMH
!+        10/05/2010
!+        V6P0
!+   CASE WITH DEPOSITION REWRITTEN, TESTS CHANGED.
!
!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 (LNHE)
!+        12/04/2011
!+        V6P0
!+   One bug corrected in case of restart, and a formula made clearer
!+   Look for 12/04/2011...
!
!history  C. VILLARET (EDF-LNHE), P.TASSI (EDF-LNHE)
!+        19/07/2011
!+        V6P1
!+  Name of variables
!+
!history  C. VILLARET (EDF-LNHE)
!+        18/01/2012
!+        V6P1
!+  Dimension 10 of AVAIL changed into NOMBLAY
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!history  R.KOPMANN (BAW)
!+        Dezember 2017
!+        V7P2
!+   Zeroed the layer thicknesses if a layer was destroyed and setting
!+   the first class to 1 and all others to zero
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
!| AVAIL          |<->| SEDIMENT FRACTION FOR EACH LAYER, CLASS, POINT
!| CONST_ALAYER   |-->| CONSTANT ACTIVE LAYER THICKNESS OR NOT
!| ELAY           |<->| ACTIVE LAYER THICKNESS FOR EACH POINT
!| ELAY0          |<->| ACTIVE LAYER THICKNESS
!| ES             |<->| LAYER THICKNESSES AS DOUBLE PRECISION
!| ESTRAT         |<->| ACTIVE STRATUM THICKNESS FOR EACH POINT
!| ESTRATNEW      |<->| ACTIVE STRATUM THICKNESS AT TIME T+DT
!| MASBAS         |-->| INTEGRAL OF TEST FUNCTIONS
!| NLAYER         |<--| NUMBER OF LAYER FOR EACH POINT
!| NLAYNEW        |<->| NUMBER OF LAYER AT TIME T+DT
!| NPOIN          |-->| NUMBER OF POINTS
!| NSICLA         |-->| NUMBER OF SIZE CLASSES FOR BED MATERIALS
!| VOLTOT         |<->| TOTAL VOLUME OF SEDIMENT IN THE BED
!| ZF             |-->| ELEVATION OF BOTTOM
!| ZFCL_W         |-->| BED EVOLUTION FOR EACH SEDIMENT CLASS
!| ZR             |-->| NON ERODABLE BED
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_sisyphe , ONLY : nomblay
      USE declarations_special
      USE interface_parallel, ONLY : p_dsum,p_isum
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      type(bief_obj),  INTENT(IN)    :: zfcl_w,zr,zf
      type(bief_obj),  INTENT(IN)    :: masbas,acladm
      INTEGER,          INTENT(IN)    :: NSICLA,NPOIN
      LOGICAL,          INTENT(IN)    :: CONST_ALAYER
      type(bief_obj),  INTENT(INOUT) :: nlayer,estrat,elay
      DOUBLE PRECISION, INTENT(INOUT) :: ELAY0
      DOUBLE PRECISION, INTENT(INOUT) :: ES(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: AVAIL(npoin,nomblay,nsicla)
      DOUBLE PRECISION, INTENT(INOUT) :: VOLTOT(nsicla),ESTRATNEW(npoin)
      INTEGER         , INTENT(INOUT) :: NLAYNEW(npoin)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,J,K,ARRET,ARRET2
      DOUBLE PRECISION EVOL,HEIGH,TEST1,TEST2,AUX
!
!-----------------------------------------------------------------------
!
!     TO CHECK FRACTIONS IN THE RANGE [-ZERO,1+ZERO]
!
!     DATA             ZERO/1.D-10/
!     IN CASE OF RESTART, THE FIRST TIME STEP IS A BIT HARD BECAUSE OF
!     SINGLE PRECISION, WITHOUT RESTART 1.D-10 IS OK
      DOUBLE PRECISION, PARAMETER :: ZERO = 1.d-7
!
!-----------------------------------------------------------------------
!
      arret=0
!
!-----------------------------------------------------------------------
!
      DO j=1,npoin
!
!       ACTIVATE BEFORE INVESTIGATING PROBLEMS IN LAYER...
!       LOOK FOR 'ACTIVATE' TO SEE OTHER LINES LIKE THIS BELOW
!
!       IF(ELAY%R(J).LT.0.D0) THEN
!         WRITE(LU,*) 'NEGATIVE ELAY IN LAYER J=',J,' ELAY=',ELAY%R(J)
!         CALL PLANTE(1)
!         STOP
!       ENDIF
!
        IF(.NOT.const_alayer) elay0 = 3.d0 * acladm%R(j)
!
        nlaynew(j) = nlayer%I(j)
!
!       TODO: QUESTION JMH, EVOLUTION HAS BEEN COMPUTED BEFORE IN ARRAY E, WHY NOT
!               EVOL=E(J) ?????
!               ELAY(J) = ES(J,1) WHY IS IT AN EXTRA ARRAY ??
!
!
        heigh = zf%R(j)-zr%R(j)
!
!       ACTIVATE BEFORE INVESTIGATING PROBLEMS IN LAYER...
!
!       IF(HEIGH.LT.0.D0) THEN
!         WRITE(LU,*) 'BAD DATA IN LAYER J=',J,' HEIGH=',HEIGH
!         CALL PLANTE(1)
!         STOP
!       ENDIF
!
!       HERE ELAY.NE.HEIGH BECAUSE ELAY IS THE ACTIVE LAYER THICKNESS
        evol  = 0.d0
        DO i=1,nsicla
          evol = evol + zfcl_w%ADR(i)%P%R(j)
        ENDDO
!
        IF(nlayer%I(j).GT.1) THEN
!
          IF(evol.GE.0.d0) THEN
!
!           DEPOSITION
!
!           NEW HEIGHT OF LAYER 2 (IT RECEIVES EVOL TO KEEP LAYER 1 CONSTANT)
            estratnew(j) = estrat%R(j) + evol
!
            DO i=1,nsicla
!             THE OLD IMPLEMENTATION CONSISTED OF FIRST
!             GIVING EVOL TO LAYER 2, WITH OLD AVAIL, THEN OF RECEIVING
!             THE DEPOSITION, BUT IF A CLASS DISAPPEARS IN LAYER 1,
!             IT IS NOT POSSIBLE TO GIVE IT FIRST TO LAYER 2, SO NEW
!             FRACTIONS MUST BE COMPUTED BEFORE GIVING EVOL TO LAYER 2
!             THEN I SEE NO DIFFERENCE BETWEEN EVOLELAY0
!             SO BOTH ARE TREATED BELOW, UNLIKE RELEASE 5.9.
!
!             1) LAYER 1 RECEIVES ZFCL_W OF CLASS I, WE COMPUTE THE
!                PROVISIONAL NEW AVAIL(J,1,I) IN AUX
              aux=(avail(j,1,i)*elay0+zfcl_w%ADR(i)%P%R(j))/
     &                        (elay0+evol)
!
!             2) LAYER 2 RECEIVES AUX*EVOL OF CLASS I (AUX MAY BE 0 HERE)
              avail(j,2,i)=(aux*evol+avail(j,2,i)*estrat%R(j))/
     &                               estratnew(j)
!
!             3) SEEN FROM LAYER 1, AUX*EVOL OF CLASS I HAS BEEN GIVEN
!                AND THE NEW LAYER THICKNESS IS ELAY0, HENCE THE NEW AVAIL
              avail(j,1,i)=( avail(j,1,i)*elay0-aux*evol
     &                          +zfcl_w%ADR(i)%P%R(j) )/elay0
!
! NOTE CV: CAN BE REPLACED BY
!              AVAIL(J,1,I)= AUX
!
!             OLD (AND WRONG) FORMULATION (IN IT -AVAIL*EVOL SHOULD BE -AUX*EVOL)
!             AVAIL(J,1,I)=( AVAIL(J,1,I)*(ELAY0-EVOL)
!    &                       +ZFCL_W%ADR(I)%P%R(J)     )/ELAY0
!
              IF(avail(j,1,i).GT.1.d0+zero.OR.
     &           avail(j,1,i).LT.-zero) THEN
                WRITE(lu,*) 'ERROR IN LAYER CASE 1'
                CALL plante(1)
                stop
              ENDIF
            ENDDO
!           NEW HEIGHT OF LAYER 1
            elay%R(j) = elay0
!
          ELSEIF(evol.GT.-elay0) THEN
! CV: I DON'T AGREE WITH THE COMMENT BELOW
!     WE'RE IN THE CASE : -ELAY0<EVOL<0 HENCE ELAY0>-EVOL>0
!     THICKNESS OF THE FIRST LAYER IS THEREFORE SUFFICIENT
!
!           EROSION GREATER THAN LAYER 1, WE HAVE TO DESTROY A STRATUM
!
            IF(-evol.GT.estrat%R(j)) THEN
!
!  CV: USUALLY, LAYER 2 IS VERY BROAD AND TWO LAYERS ARE IN GENERAL SUFFICIENT
!      HERE LAYER 2 IS DESTROYED
!
!             USUAL CASE (NOTE JMH : WHY NOT .GE.2 ?
!
              IF(nlayer%I(j).GT.2) THEN
!
                DO i=1,nsicla
                  avail(j,1,i) = (  avail(j,1,i)*elay0
     &                            + zfcl_w%ADR(i)%P%R(j)
     &                            + avail(j,2,i)*estrat%R(j)
     &                            - avail(j,3,i)*(evol+estrat%R(j))
     &                           )/ elay0
                  IF(avail(j,1,i).GT.1.d0+zero.OR.
     &               avail(j,1,i).LT.-zero) THEN
                    WRITE(lu,*) 'ERROR IN LAYER CASE 2'
                    CALL plante(1)
                    stop
                  ENDIF
                  avail(j,2,i) = avail(j,3,i)
                  DO k=3,min(nomblay-1,nlayer%I(j))
                    avail(j,k,i) = avail(j,k+1,i)
                  ENDDO
                ENDDO
                elay%R(j) = elay0
                nlaynew(j) = nlayer%I(j) - 1
                estratnew(j) = estrat%R(j) + evol + es(j,3)
                DO k=3,min(nomblay-1,nlayer%I(j))
                  es(j,k) = es(j,k+1)
                ENDDO
!               last layer must be zeroed, AVAIL: first class = 1, others = 0
                es(j,nlayer%I(j)) = 0.d0
                avail(j,nlayer%I(j),1) = 1.d0
                DO i=2,nsicla
                  avail(j,nlayer%I(j),i) = 0.d0
                ENDDO


!
!             HERE NLAYER.GT.1 AND NOT .GT.2, SO 2 !
!
              ELSE
                DO i=1,nsicla
                  IF(heigh.GT.0.d0) THEN
                    avail(j,1,i) = (  avail(j,1,i)*elay0
     &                              + zfcl_w%ADR(i)%P%R(j)
     &                              + estrat%R(j)*avail(j,2,i)
     &                              )/(elay0+evol+estrat%R(j))
!                   MODIF JMH 12/04/2011, WITH 2 LAYERS, LAYER 1
!                   HAS ELAY%R(J)=ELAY0, SO WHY ELAY IN DENOMINATOR
!                   AND ELAY0 IN NUMERATOR ? (VICIOUS!!)
!    &                              )/(ELAY%R(J)+EVOL+ESTRAT%R(J))
                    IF(avail(j,1,i).GT.1.d0+zero.OR.
     &                 avail(j,1,i).LT.-zero) THEN
                      WRITE(lu,*) 'J=',j,' NLAYER%I(J)=',nlayer%I(j)
                      WRITE(lu,*) 'AVAIL(J,1,I)=',avail(j,1,i)
                      WRITE(lu,*) 'AVAIL(J,2,I)=',avail(j,2,i)
                      WRITE(lu,*) 'ZFCL=',zfcl_w%ADR(i)%P%R(j)
                      WRITE(lu,*) 'HEIGH=',heigh,' ELAY0=',elay0
                      WRITE(lu,*) 'ESTRAT%R(J)=',estrat%R(j)
                      WRITE(lu,*) 'ELAY%R(J)=',elay%R(j)
                      WRITE(lu,*) 'EVOL=',evol
                      WRITE(lu,*) 'ERROR IN LAYER CASE 3'
                      CALL plante(1)
                      stop
                    ENDIF
                  ELSE
                    avail(j,1,i) = 0.d0
! Fraction 1 should be 1, if all the others are zero
                    avail(j,1,1) = 1.d0
                  ENDIF
                  avail(j,2,i) = 0.d0
                ENDDO
! Fraction 1 should be 1, if all the others are zero
                avail(j,2,1) = 1.d0
                nlaynew(j) = nlayer%I(j) - 1
                elay%R(j) = heigh
                estratnew(j) = 0.d0
!           a layer was just canceled, the layer thickness of the canceled layer
!           should be set to zero for graphic printout and restart
                es(j,nlaynew(j)+1) = 0.d0
              ENDIF
!
!           ONLY LAYER 1 ERODED
!
            ELSE
              DO i=1,nsicla
                avail(j,1,i) = (  avail(j,1,i) * elay0
     &                          + zfcl_w%ADR(i)%P%R(j)
     &                          - evol*avail(j,2,i)    )/elay0
                IF(avail(j,1,i).GT.1.d0+zero.OR.
     &             avail(j,1,i).LT.-zero) THEN
                  WRITE(lu,*) 'ERROR IN LAYER CASE 4'
                  CALL plante(1)
                  stop
                ENDIF
              ENDDO
              elay%R(j) = elay0
              estratnew(j) = estrat%R(j) + evol
            ENDIF
!
          ELSE
!
!           STOPS IF EROSION IS GREATER THAN
!           THICKNESS OF THE ACTIVE LAYER!
!
            WRITE(lu,*) 'TOO MUCH EROSION AT POINT J=',j
            WRITE(lu,*) 'DECREASE DT OR INCREASE ELAY0'
            WRITE(lu,*) 'EVOL=', evol, 'ELAY0=',elay0
            CALL plante(1)
            stop
!
!         END OF TESTS ON EVOL
!
          ENDIF
!
! THERE WAS ONLY ONE LAYER
! ------------------------
!
        ELSE
!
!         IT IS NOW BIG ENOUGH TO MAKE TWO LAYERS
!
          IF(heigh.GT.elay0) THEN
            nlaynew(j) = 2
            estratnew(j) = heigh - elay0
            elay%R(j) = elay0
            DO i=1,nsicla
              avail(j,2,i) = avail(j,1,i)
              avail(j,1,i) = (avail(j,1,i) * (elay0-evol)
     &                     + zfcl_w%ADR(i)%P%R(j) )/elay0
              IF(avail(j,1,i).GT.1.d0+zero.OR.
     &          avail(j,1,i).LT.-zero) THEN
                WRITE(lu,*) 'ERROR IN LAYER CASE 5'
                CALL plante(1)
                stop
              ENDIF
            ENDDO
!
! IF THERE REMAINS ONLY ONE LAYER
! -------------------------------
!
          ELSE
!           NOTE JMH: THE TRICKIEST PART...
!           THE PROBLEM OF 0/0 CREATED BY THE CHOICE OF AVAIL
!           AS MAIN VARIABLE...
            IF(elay%R(j)+evol.GT.1.d-15) THEN
              DO i=1,nsicla
!               AUX=AVAIL(J,1,I)
                avail(j,1,i) = (avail(j,1,i)*elay%R(j)+
     &                          zfcl_w%ADR(i)%P%R(j))
     &                          / (elay%R(j)+evol)
!               IF(AVAIL(J,1,I).GT.1.D0+ZERO.OR.
!    &            AVAIL(J,1,I).LT.-ZERO) THEN
!                 WRITE(LU,*) 'ERROR IN LAYER CASE 6'
!                 WRITE(LU,*) 'INITIAL AVAIL=',AUX
!                 WRITE(LU,*) 'J=',J,' CLASS ',I
!                 WRITE(LU,*) 'EVOL=',EVOL,' ELAY=',ELAY%R(J)
!                 WRITE(LU,*) 'EVOL+ELAY=',EVOL+ELAY%R(J)
!                 WRITE(LU,*) 'ZFCL=',ZFCL_W%ADR(I)%P%R(J)
!                 WRITE(LU,*) 'DENOMINATOR=',ELAY%R(J)+EVOL
!                 WRITE(LU,*) 'NUMERATOR=',AUX*ELAY%R(J)+
!    &                                     ZFCL_W%ADR(I)%P%R(J)
!               ENDIF
                avail(j,2,i) = 0.d0
              ENDDO
              IF(elay%R(j)+evol.LT.1.d-5) THEN
!               PLAYING WITH ZEROES, RISK OF SUM NOT EQUAL TO 1
!               ONLY BECAUSE OF TRUNCATION ERRORS, WE NORMALISE
                test1=0.d0
                DO i=1,nsicla
                  avail(j,1,i)=max(0.d0,min(1.d0,avail(j,1,i)))
                  test1=test1+avail(j,1,i)
                ENDDO
                IF((test1-1.d0)**2.GT.zero) THEN
                  DO i=1,nsicla
                    avail(j,1,i)=avail(j,1,i)/max(test1,1.d-21)
                  ENDDO
                ENDIF
              ENDIF
            ELSE
!             ACLADM BEING RECOMPUTED WITH AVAIL
              avail(j,1,1)=1.d0
              avail(j,2,1)=1.d0
              DO i=2,nsicla
                avail(j,1,i) = 0.d0
                avail(j,2,i) = 0.d0
              ENDDO
            ENDIF
            elay%R(j) = heigh
            estratnew(j) = 0.d0
            nlaynew(j) = 1
!RK
            DO i=2,nomblay
              es(j,i) = 0.d0
            END DO
!
          ENDIF
        ENDIF
!
      nlayer%I(j) = nlaynew(j)
      estrat%R(j) = estratnew(j)
      es(j,1) = elay%R(j)
!     IN CASE OF RESTART, ES(J,2)
!     WILL BE STORED IN A FILE AND LOOKED AT TO COUNT THE
!     NUMBER OF LAYERS
      es(j,2) = estrat%R(j)
!
      test1 = 0.d0
      test2 = 0.d0
!
      DO i=1,nsicla
        DO k=1,nlayer%I(j)
!         CHECKS THAT AVAIL IS IN THE RANGE (-ZERO,1+ZERO)
          IF(avail(j,k,i).GT.1.d0+zero.OR.avail(j,k,i).LT.-zero) THEN
            WRITE(lu,*) 'ERROR ON FRACTIONS'
            WRITE(lu,*) 'LAYER ',k,' CLASS ',i,' POINT ',j
            IF(avail(j,k,i).LT.0.d0) THEN
              WRITE(lu,*) 'AVAIL=' ,avail(j,k,i)
            ELSE
              WRITE(lu,*) 'AVAIL-1=' ,avail(j,k,i)-1.d0
            ENDIF
            WRITE(lu,*) 'ZFCL=',zfcl_w%ADR(i)%P%R(j)
            WRITE(lu,*) 'EVOL=',evol,' ELAY=',elay%R(j)
            arret=1
          ELSE
!           ONCE CHECKED THAT WE HAVE ONLY TRUNCATION ERRORS, CLIPS
            avail(j,1,i)=max(avail(j,1,i),0.d0)
            avail(j,1,i)=min(avail(j,1,i),1.d0)
          ENDIF
        ENDDO
        test1 = test1 + avail(j,1,i)
        test2 = test2 + avail(j,2,i)
      ENDDO
!
!     CHECKS THAT SUM OF AVAIL IS 1 FOR FIRST 2 LAYERS
!
      IF(test1.GT.zero.AND.(test1-1.d0)**2>zero) THEN
        WRITE(lu,*) ' PROBLEM IN LAYER: J,TEST1',j,test1
        WRITE(lu,*) ' IN LAYER 1 SUM OF FRACTIONS NOT 1'
        arret=1
      ENDIF
      IF(test2.GT.zero.AND.(test2-1.d0)**2>zero) THEN
        WRITE(lu,*) ' PROBLEM IN LAYER: J,TEST2',j,test2
        WRITE(lu,*) ' IN LAYER 2 SUM OF FRACTIONS IS NOT 1'
        arret=1
      ENDIF
!
!     END OF LOOP ON ALL POINTS
!
      ENDDO
!
!     COMPUTES THE TOTAL VOLUME OF SEDIMENTS IN THE DOMAIN
!
      DO i = 1, nsicla
        voltot(i) = 0.d0
      ENDDO
      DO i=1,nsicla
        DO j=1,npoin
          DO k=1,nlayer%I(j)
            voltot(i) = voltot(i) + es(j,k)*avail(j,k,i)*masbas%R(j)
          ENDDO
        ENDDO
      ENDDO
!
      IF(ncsize.GT.1) THEN
        DO i=1,nsicla
          voltot(i) = p_dsum(voltot(i))
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
!     CLEAN STOP FOR ALL PROCESSORS IF PROBLEM
!
      arret2=arret
      IF(ncsize.GT.1) arret2=p_isum(arret)
      IF(arret2.GT.0) THEN
        WRITE(lu,*) 'STOP AFTER AN ERROR IN LAYER'
        IF(arret.EQ.0) THEN
          WRITE(lu,*) 'IN ',arret2,' PROCESSOR(S)'
        ENDIF
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END