init_avai.f

Go to the documentation of this file.

!                   ********************
                    SUBROUTINE init_avai
!                   ********************
!
!
!***********************************************************************
! SISYPHE   V6P2                                   21/07/2011
!***********************************************************************
!
!brief    INITIAL FRACTION DISTRIBUTION AND LAYER THICKNESS.
!
!note     PHILOSOPHY: INIT_COMPO DEFINES THE STRATIFICATION CORRESPONDING
!+                TO THE BOTTOM INITIAL COMPOSITION; NCOUCHES CORRESPONDS
!+                TO THE NUMBER OF REAL INITIAL LAYERS.
!note         INIT_AVAI CORRECTS AND SUPPLEMENTS THIS STRATIFICATION
!+                IF PROBLEMS, AND ADDS THE ACTIVE LAYER; NLAYER CORRESPONDS
!+                TO THE NUMBER OF LAYERS USED IN THE COMPUTATION.
!
!history  BUI MINH DUC
!+        2002
!+
!+   INITIAL FRACTION DISTRIBUTION FOR NON-UNIFORM BED MATERIALS
!
!history  MATTHIEU GONZALES DE LINARES
!+        2002/2003
!+
!+
!
!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.RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!!history  R.KOPMANN (BAW)
!+        Januar 2018
!+        V7P3
!+   Checks for bad layer thicknesses and fractions
!
!history  R.KOPMANN (BAW)
!+        15/02/2019
!+        V7P2
!+   Adding initial volume for mass balance
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_telemac
      USE declarations_sisyphe
      USE interface_sisyphe, ex_init_avai=> init_avai
!
      USE declarations_special
      USE interface_parallel, ONLY : p_dsum
      IMPLICIT NONE
!
!
!-----------------------------------------------------------------------
!
      INTEGER I,J,K,NMAXI
!
      DOUBLE PRECISION HAUTSED,TEST1
!
!-----------------------------------------------------------------------
!
      nmaxi = 0
!
!     THE INITIAL NUMBER OF LAYERS, THEIR THICKNESS AND THEIR COMPOSITION
!     ARE SET BY THE USER
!
!     NOTE: WHEN COMPUTATION CONTINUED INIT_COMPO MUST NOT
!           CHANGE ES AND AVAIL
!
      IF(debu) THEN
!       TODO: TENTATIVE VALUE, THIS IS TO CHECK
        DO j=1,npoin
          i=1
          IF(nomblay.GE.2) THEN
            DO k=2,nomblay
              IF(es(j,k).GT.0.d0) i = i + 1
            ENDDO
          ENDIF
          nlayer%I(j)=i
!         CHECKING ALL LAYERS AND CORRECTING AVAIL
!         DUE TO POSSIBLE SHIFT OF SINGLE PRECISION STORAGE
          DO k=1,nlayer%I(j)
            test1=0.d0
            DO i=1,nsicla
              test1=test1+avail(j,k,i)
            ENDDO
            IF(test1.GT.1.d-10.AND.(test1-1.d0)**2.GT.1.d-10) THEN
              DO i=1,nsicla
                avail(j,k,i)=avail(j,k,i)/test1
              ENDDO
            ENDIF
          ENDDO
        ENDDO
      ELSE
!
        CALL init_compo(it1%I)
!
        DO j=1,npoin
!
!       NOMBLAY IS THE MAXIMUM NUMBER OF LAYERS ALLOWED
!
        nlayer%I(j) = it1%I(j)
        IF(nlayer%I(j).GT.nomblay) THEN
          WRITE(lu,1815) nomblay
          CALL plante(1)
          stop
        ENDIF
!
!       THE HEIGHT OF SEDIMENT (SUM OF ES) MUST NOT BE MORE THAN ZF-ZR
!       IF SO, THE HEIGHT OF THE LAST LAYER IS REDUCED
!       IF THERE ARE LAYERS UNDER ZR, THEY ARE NOT TAKEN INTO ACCOUNT
        hautsed = 0.d0
        DO k=1,it1%I(j)
          IF(hautsed + es(j,k) .GE. zf%R(j) - zr%R(j)) THEN
            es(j,k) = zf%R(j) - zr%R(j) -  hautsed
            nlayer%I(j) = k
            hautsed = hautsed + es(j,k)
            EXIT
          ENDIF
          hautsed = hautsed + es(j,k)
        ENDDO
!
!       OTHER LAYERS SET TO 0.D0
!
        IF(nlayer%I(j).LT.nomblay) THEN
          DO k=nlayer%I(j)+1,nomblay
            es(j,k) = 0.d0
          ENDDO
        ENDIF
!
!       THE THICKNESS OF THE LAST LAYER IS ENLARGED SO THAT
!       THE HEIGHT OF SEDIMENT (SUM OF ES) IS EQUAL TO ZF-ZR
!
        IF(hautsed.LT.zf%R(j)-zr%R(j)) THEN
          es(j,nlayer%I(j))=es(j,nlayer%I(j))+zf%R(j)-zr%R(j)-hautsed
        ENDIF
!
        IF(nlayer%I(j).GT.1) THEN
!         IT IS ASSUMED THAT ELAY0 IS SMALLER THAN THE FIRST STRATUM
!         NEED TO ADD THE CASE WHERE ELAY0 IS LARGER
          IF(elay0.GT.es(j,1)) THEN
            WRITE(lu,*) ' ACTIVE LAYER TOO BIG/STRATIFICATION'
            CALL plante(1)
            stop
          ENDIF
        ENDIF
!
!       THE FIRST STRATUM IS SEPARATED BETWEEN ACTIVE LAYER + ACTIVE STRATUM
!
        IF(elay0.LT.es(j,1)) THEN
          nlayer%I(j) = nlayer%I(j) + 1
          IF(nlayer%I(j).GT.nomblay) THEN
            WRITE(lu,1815) nomblay
            CALL plante(1)
            stop
          ENDIF
!         INDICES FOR ES AND AVAIL NEED TO BE OFFSET
          IF(nlayer%I(j).GT.2) THEN
            DO k=nlayer%I(j),3,-1
              es(j,k) = es(j,k-1)
            ENDDO
          ENDIF
          es(j,2) = es(j,1) - elay0
          es(j,1) = elay0
          DO i=1,nsicla
            DO k=nlayer%I(j),2,-1
              avail(j,k,i) = avail(j,k-1,i)
            ENDDO
          ENDDO
        ENDIF
!
        ENDDO !J
!
      ENDIF
!
      nmaxi=0
      DO j=1,npoin
!
!       CHECKS FOR BAD LAYER THICKNESSES
        DO k=1,nomblay
          IF(es(j,k).LT.0.d0) THEN
            WRITE(lu,*)'BAD LAYER THICKNESS',j,k ,' :',es(j,k)
            CALL plante(1)
            stop
          ENDIF
!       CHECKS FOR BAD AVAI
          DO i=1,nsicla
            IF(avail(j,k,i).LT.0.d0.OR.avail(j,k,i).GT.1.d0) THEN
              WRITE(lu,*)'BAD FRACTIONS',j,k,i ,' :',avail(j,k,i)
              CALL plante(1)
              stop
            ENDIF
          END DO
        END DO
!
        elay%R(j) = es(j,1)
        IF(nlayer%I(j).GT.1) THEN
          estrat%R(j) = es(j,2)
        ENDIF
!
!       UNUSED AVAIL ARE FILLED WITH ZEROS (IS IT USEFUL ???)
!
        IF(nlayer%I(j).LT.nomblay) THEN
          DO i = 1, nsicla
            DO k = nlayer%I(j)+1,nomblay
              avail(j,k,i) = 0.d0
            ENDDO
          ENDDO
        ENDIF
        IF(nlayer%I(j).GT.nmaxi) nmaxi = nlayer%I(j)
      ENDDO
!
!     COMPUTES THE TOTAL VOLUME OF SEDIMENTS IN EACH CLASS
!
      DO i=1,nsicla
        voltot(i)=0.d0
        DO j=1,npoin
          DO k=1,nlayer%I(j)
            voltot(i) = voltot(i) + es(j,k)*avail(j,k,i)*volu2d%R(j)
          ENDDO
        ENDDO
      ENDDO
!
      IF(ncsize.GT.1) THEN
        DO i=1,nsicla
          voltot(i) = p_dsum(voltot(i))
        ENDDO
      ENDIF
      volini = voltot
!
      WRITE(lu,*) 'MAXIMUM INITIAL NUMBER OF LAYERS :',nmaxi
      DO i=1,nsicla
        WRITE(lu,*)'TOTAL VOLUME OF CLASS ',i ,' :',voltot(i)
      ENDDO
!
!-----------------------------------------------------------------------
!
1815  FORMAT(1x,'THERE ARE MORE THAN ',1i3,' LAYERS IN STRATIFICATION')
!
!-----------------------------------------------------------------------
!
      RETURN
      END